3.1 Elliptic PDEs and operators

Inverse problems studied by Factas involve systems governed by an equation of the form $ℒ\varphi =\Psi \left(m\right)$, where $ℒ$ is an elliptic partial differential operator and $\Psi \left(m\right)$ a source term depending on some unknown quantity $m$. The data consist of incomplete measurements of the potential $\varphi$ or its gradient (the field) in a portion of space, away from the support of the source.

3.2 Inverse source problems

Given an elliptic PDE of the form $ℒ\varphi =\Psi \left(m\right)$ as in Section 3.1, the corresponding inverse source problem consists in recovering the quantity $m$ from the data, that typically consist of measurements of the potential $\varphi$ or the field away from the support of the source. Usually the support of $m$ is assumed to be compact, and a super-set thereof is specified. This kind of issues may thus be seen as parameterized inverse potential problems (parameterized by $m$, that is). They arise naturally in non-destructive testing, medical imaging, paleo-magnetism, gravimetry and geosciences, as well as in inverse scattering, see Section 4. The second and third application domains pertain to static electromagnetism, a framework in which source terms typically occur in divergence form; that is, $\Psi \left(m\right)=\mathrm{div}m$ where $m$ is a ${ℝ}^3$-valued field or distribution. As a rule, the operator $ℒ$ is then of the form $ℒ\varphi =\mathrm{div}\left(\Sigma \nabla \varphi \right)$, where $\Sigma$ is valued in positive definite matrices satisfying fixed ellipticity bounds and relates to the electromagnetic characteristics of the medium. In this case, the problem amounts to recover a vector field (namely: $m/\det \left(\Sigma \right)$) knowing a super-set $S$ of its support and the gradient summand of its Helmholtz decomposition outside of $S$, for the Riemannian metric with tensor $\left(\det \Sigma \right){\Sigma }^{-1}$. The simplest case, of course, occurs in the Euclidean setting, where $ℒ$ is the ordinary Laplacian, which corresponds to homogeneous media.

3.3 Rational approximation, behavior of poles

Rational approximation to holomorphic functions of one complex variable is a long standing chapter of classical analysis, with notable applications to number theory, spectral theory and numerical analysis. Over the last decades, it has become a cornerstone of modeling in Systems Engineering, and it can also be construed as a technique to regularize inverse source problems in the plane, where the degree is the regularizing parameter. Indeed, by partial fraction expansion, a rational function can be viewed as the complex derivative of a discrete logarithmic potential with as many masses as the degree (assuming that the poles are simple); that is, if $f$ is rational of degree $N$, then $\overline{\partial }f={\sum }_{j=1}^N{a}_j{\delta }_{z_j}$ where the $z_j$ are its poles and ${\delta }_{z_j}$ is a Dirac unit mass at $z_j$. Moreover, a holomorphic function is the complex derivative of the logarithmic potential of its own values on a curve encompassing the domain of analyticity (this is the Cauchy formula); hence, rational approximation aims at representing as well as possible a logarithmic potential by a discrete potential with prescribed number of masses, in the sense that their derivatives should be close (a Sobolev-type approximation).

Predecessors of Factas (the Apics and Miaou project teams) have designed a dedicated steepest-descent algorithm for quadratic approximation criteria whose convergence to a local minimum is guaranteed. This gradient algorithm may either be initialized by a preliminary approximation method, or recursively proceed with respect to the degree $N$ of the approximant, on a compactification of the parameter space 35, as can be done with the RARL2 software (see Section 3.5.3). It has proved to be effective in applications carried out by the team (see for instance 12 for the identification of micro-wave filters, and Sections 4.1, 4.3).

However, finding best rational approximants of prescribed degree to a specific function, say in the uniform norm on a given set, seems out of reach except in rare, particular cases. Instead, constructive rational approximation has focused on estimating optimal convergence rates and deriving approximation schemes coming close to meet them, or studying computationally appealing approximants like Padé interpolants and their variants. Two main issues are then the effective computation of optimal or near optimal approximants of given degree, and the connection between the singularities of the approximant (the poles) and those of the approximated function. Factas has contributed to both.

As regards near-optimal approximants, their design requires a knowledge of optimal rates in the situation at hand. In recent years, we were active determining lower bounds on that rate, a piece of information which is crucial but difficult to obtain. Our methods are topological in nature (Ljustenik-Schnirelman theory, genus of compact symmetric sets), like most techniques in the area, and in collaboration with Q. Tao from the University of Macao, we devised algorithms to compute lower bounds in best $L^2$ approximation by stable rational function of given degree on the unit circle, which is first of this kind and sometimes quite precise, see 4. Research in this direction is still going on, in particular on best $L^2$ approximation of uni-modular functions that cannot be approximated in sup norm, except when they are rational of admissible degree, in which case they are, of course, their own best approximant.

We refer here to the behavior of best rational approximants of given degree, in the $L^{\infty }$-sense on a compact subset of the domain of analyticity, to complex analytic functions $f$ that can be continued analytically (possibly in a multi-valued manner) except perhaps over a set of logarithmic capacity zero in the plane. When the continuation of $f$ has finitely many branches; that is, if the Riemann surface to which this continuation extends analytically in a single-valued manner (except perhaps on a polar subset thereof) is compact, then the behavior as $N$ goes large of rational approximants of degree $N$ whose $N$-th root error is asymptotically smallest possible (in particular the asymptotic behavior of best approximants) has recently been elucidated rather completely in this joint research effort involving Factas.

Another classical technique to approximate –more accurately: extrapolate– a function, given a set of pointwise values, is to compute a rational interpolant of minimal degree to match the values. This method, known as Padé (or multi-point Padé) approximation has been intensively studied for decades 32 but fails to produce pointwise convergence, even if the data are analytic. The best it can give in general, at least to functions whose singular set has capacity zero, is convergence in capacity which does not prevent poles of the approximant from wandering about the domain of analyticity of the approximated function, but does imply that each pole of the approximated function attracts a pole of the approximant 71. This phenomenon is well-known in numerical analysis, and leads Physicists and Engineers to distinguish between “mathematical” and “physical” poles. A modification of the multi-point Padé technique, in which the degree is kept much smaller than the number of data and approximate interpolation is performed in the least-square sense, has become especially popular over the last decade under the name vector fitting; it teams up with a barycentric representation of rational functions satisfying prescribed interpolation conditions, known as AAA (for Anderson-Antoulas Adaptive) scheme. Though the behavior of this least square substitute to Padé approximation, defined by Equation (1) in Section 4.4, resembles the one of multi-point Padé approximants from a numerical viewpoint, there has been apparently no convergence result for such approximate interpolants so far. Motivated by the outcome of numerical schemes developed by our partners to recover resonance frequencies of conductors under electromagnetic inverse scattering, the PhD thesis 31 of P. Asensio started investigating the behavior of such least-square rational approximants to functions with polar singular set, and dwelling on this work, we were recently able to show convergence in capacity thereof.

Regarding complex rational approximation as a means to tackle inverse source problems in the plane makes for a unifying point of view on various deconvolution techniques, from system identification and time series analysis to frequency-wise inverse scattering and non-destructive testing. But still more interestingly perhaps, it is suggestive of similar approaches to problems in higher dimension, where holomorphic functions generalize to harmonic gradients and rational functions to finite linear combinations of dipoles, see Section 3.1.2. This line of research is only starting, but seems to offer new avenues in connection with applications.

3.4 Asymptotic analysis

Asymptotic analysis deals with understanding behavior or explicit construction of the solution when a parameter entering a problem is either small or large. Factas has been involved in applications of asymptotic analysis in different contexts including both formal constructions and their rigorous justifications.

One type of asymptotic analysis for dynamical problems is the large-time behavior analysis. A rather classical issue here is that of limiting amplitude principle for wave equation. This principle states that the solution of the time-dependent wave equation with a periodic-in-time monochromatic source term $f{e}^{i\omega t}$ necessarily stabilizes for large times $t$ to the solution of the corresponding Helmholtz equation: $\mathrm{div}\left(\alpha \nabla \varphi \right)+{\omega }^2\beta \varphi =-f$ for suitable known functions $\alpha$, $\beta$ and $f$ depending on space variables. Revisiting this topic is motivated by recent popularity of numerical time-domain approaches (e.g., 28, 29, 64, 86) to elliptic PDE problems through efficient solution of auxiliary time-dependent equations where, for example, computational effort needed to be concentrated only on wave front neighborhood which is small for high-frequencies, a notoriously difficult regime for numerical solution of Helmholtz problems. With this respect, not only the fact of the time convergence in the limiting amplitude principle is important but also quantification of its rate. Dmitry Ponomarev was involved in the work 13 that deals with the limiting amplitude principle for a non-homogeneous medium wave-equation in different dimensions (in the dimension 1, a slight modification of the classical limiting amplitude principle was proposed). The analysis approach relies on proving that the solution decays to zero for a source-free problem with localized initial data and radiation boundary condition, an interesting problem of asymptotic analysis on its own, even in dimension 1 30. In 80, convergence to the periodic motion was also shown in a totally different context of one model of mechanical sliding contact problem with wear 20.

Previously, in a nonlinear context, rigorous asymptotic analysis 77 was instrumental to justify a parabolic model of pulse propagation in photo-polymers by comparing solutions of that model with those of the original Maxwell's system.

In the context of inverse problems, asymptotic analysis is useful when applied to the magnitude of the regularization parameter. When the latter tends to its limiting value, a solution of the regularized problem with ideal (noiseless) input data should tend to the exact solution. In presence of noise, it is important to relate this convergence rate to the value of the problem's constraint in the asymptotic regime of the regularization parameter.

The works 46, 79 on convolution integral equations on large domains, mentioned in Section 3.1.3, are an example of constructive asymptotic analysis. Here, one of the difficulty comes from a singular perturbation. Indeed, in the asymptotic limit of infinite size of the region, the spectral problem solution cease to exist since the integral operator loses the compactness property.

In the context of net magnetization reconstruction in the inverse magnetization problem (Sections 3.2 and 4.2), situation when the measurement area size is large leads to a different kind of application of constructive asymptotic analysis 36, 42, 81. Here, explicit constructions of the solution estimates are performed to different asymptotic orders with respect to the measurement region size. The higher-order estimates can give good accuracy already for relatively small value of the measurement region but are much more unstable with respect to the perturbation of the measured data. This problem also exhibits another interesting asymptotic phenomenon which is somewhat similar to the “boundary layer” common for boundary-value problem for differential equations with small or large parameters. In particular, the solution (for tangential components of net moment) is composed of a global leading-order quantity (where formal passage to the asymptotic limit can be performed) and a correction term which is localized in a region that shrinks in the asymptotic limit.

3.5 Software tools of the team

In addition to the above-mentioned research activities, Factas develops and maintains a number of long-term software tools that either implement and illustrate effectiveness of the algorithms theoretically developed by the team or serve as tools to help further research by team members. We briefly present the most important of them, which have been developed over the past few years.

4.1 Some inverse problems for cerebral imaging

Solving over-determined Cauchy problems for the Laplace equation on a spherical layer (in 3-D) in order to extrapolate incomplete data (see Section 3.1.1) is an ingredient of the team's approach to inverse source problems in electro-encephalography (EEG), see 9. The model comes from Maxwell's equation in the quasi-static regime, whence $\mathrm{div}\left(\sigma \nabla \varphi \right)=\mathrm{div}m$ in a ball $\Omega$ which is the union of nested layers ${\Omega }_i$ (brain, skull, scalp for $i=0,1,2$), where the singularities lie in ${\Omega }_0$ (i.e., the current sources $m$ are supported in $S\subset {\Omega }_0$, with the notation of Section 3), see Figure51. The inverse EEG source problem consists in recovering $m$ from pointwise values of the electric potential $\varphi$ measured on the scalp ${\Gamma }_2$ (together with the assumption that the normal current vanishes there). It first involves transmitting the data from the scalp ${\Gamma }_2$ down to the cortex ${\Gamma }_0$. Whenever the ${\Omega }_i$ are of different constant conductivities, $\varphi$ satisfies Laplace equation in the outermost shells ${\Omega }_2$ and ${\Omega }_1$ and this “cortical mapping” step is performed using integral representations of $\varphi$ on the spheres ${\Gamma }_i$ (obtained through convolution by the Poisson kernels of the balls and using Green formula, a specific use of boundary element methods), followed by expansions on spherical harmonic bases and Tikhonov regularization.

Assuming $m$ to be a linear combination of Dirac masses (dipolar sources), it turns out (by convolution with the fundamental solution of 3-D Laplace equation) that traces of $\varphi$ on 2-D cross sections of ${\Gamma }_0$ coincide with functions with branched singularities in the slicing plane 7, 43. These singularities are related to the actual location of the sources. Hence we are back to the 2-D framework of Section 3.3, and recovering these singularities can be performed via best rational approximation. The goal is to produce a fast and sufficiently accurate initial guess on the number and location of the sources in order to run heavier descent algorithms on the direct problem, which are more precise but computationally costly and often fail to converge if not properly initialized. Such a localization process can add a geometric, valuable piece of information to the standard temporal analysis of EEG signal records. It appears that, in the rational approximation step, multiple poles possess a nice behavior with respect to branched singularities. This is due to the very physical assumptions on the model from dipolar current sources, for both EEG data and MEG (magneto-encephalography) data that correspond to measurements of the magnetic field, as well as for (magnetic) field data produced by magnetic dipolar sources within rocks (see Section 4.2). Though numerically observed in 9, there is no mathematical justification so far why branched singularities generate such strong accumulation of the poles of the approximants. This intriguing property, however, is definitely helping source recovery and will be the topic of further study. It is used in order to automatically estimate the “most plausible” number of sources (numerically: up to 3, at the moment). In this connection, a software FindSources3D (FS3D, see Section 3.5.4) dedicated to pointwise source estimation in EEG–MEG has been developed. We also studied the uniqueness of a critical point of the quadratic criterion in the EEG source problem in ${\Omega }_0$ for a single dipole situation (see 15), an important issue for the use of descent algorithms.

Together with M. Darbas (LAGA, Université Sorbonne Paris Nord) and P.-H. Tournier (JLL laboratory, Sorbonne Université), we recently considered the EEG inverse problem with a variable conductivity in the intermediate skull layer ${\Omega }_1$, in order to model hard / spongy bones, especially for neonates. The related transmission step is then performed using a mixed variational regularization and finite elements on tetrahedral meshes, and the coupling with FS3D for dipolar source estimation furnishes promising results 56.

Other approaches have been studied for EEG, MEG and “Stereo” EEG (SEEG), where the potential is measured by deep electrodes and sensors within the brain as in the scheme of Figure 1, and for more realistic geometries of the head.

Show image description

Three nested ovoids illustrating the three layers of the head (denoted ${\Omega }_i$ with $i=0,1,2$, 0 being the index of the inner layer). The conductivity of the $i$-th layer is denoted with ${\sigma }_i$ and its boundary is denoted with ${\Gamma }_i$. Two red straight lines illustrate how deep electrodes take measurements in the inner layer. Little segments perpendicular to each electrode illustrate the different sensors that lie along the electrode. A blue closed curve, denoted $S$, in the inner layer illustrates the surface where the sources lie, and little arrows show that their orientation is normal to the curve.

Assuming that the current source term $m$ is a ${ℝ}^3$-valued vector field (measure, or distribution) supported on a surface $S\subset {\Omega }_0$ (the gray / white matters interface within the brain) and normally oriented to $S$, they consist in regularizing the inverse source problem by a total variation (TV) constraint - to favor sparsity - on $m$, added to the quadratic data approximation criterion (see Section 3.2). This is similar to the path that is taken for inverse magnetization problems (see Section  4.2). The approach follows that of 8 and is implemented through algorithms whose convergence properties were studied in 14, 75. We are now able to handle MEG, EEG, SEEG modalities, simultaneously or not. The simultaneous handling of the different modalities is made more straightforward by coupling the source localization problem and the inverse transmission problem. Tests on synthetic data provided good quality results, though they are quite numerically costly to obtain. This opens up the possibility to consider sources that may exhibit properties usually associated with distributions rather than functions.

4.2 Inverse magnetization issues for planar and volumetric samples

Among other things, geoscientists are interested in understanding the magnetic characteristics of ancient rocks. Indeed, ferro-magnetic particles in a rock carry a magnetization that has been acquired when the rock was hot, under the influence of the magnetic field that was ambient at that time. For an igneous rock for instance, and if no subsequent event has heated it up, this corresponds to the time when the rock was formed. If the rock can be dated, recovering its magnetization hence provides valuable information about the history of the magnetic field. This gives elements for better understanding to key questions such as: what was the magnetic field of the sun when the solar system was at the proto-planetary phase? when did the magnetic dynamo of Mars stop? when did the magnetic dynamo of Earth start?

The magnetization of a rock is not directly measurable. However, it produces a tiny magnetic field, which can be measured if the sample is isolated from other sources of magnetic field. A category of instruments of particular importance with that respect are the magnetic microscopes. They are used to measure the field produced by a fairly small sample: either a simple grain, or a wider sample that has been first prepared by gluing it on some support and polishing it until getting only a thin slab. The microscope operates at some distance above the sample and measures the magnetic field. The typical experimental set up is represented on Figure 2.

Show image description

Schematic view of the experimental setup : the horizontal plane on which the rock sample lies (the sample being a parallelepiped with height $r$ and a square basis of half-size $s$). Its basis is at height 0, and the center of the square basis is at the origin of the system of coordinates). Above, the measurements are performed on a horizontal plane at height $z$. Here, for instance, it is a square, whose half-size is $R$ and which is horizontally centered.

The magnetization $m$ is modeled as a ${ℝ}^3$-valued vector field defined on the rock sample which is assumed to be a subset of $S$. One advantage of the magnetic microscopes is that they operate close to the sample, i.e., the height $z$ of measurement is small compared to the horizontal characteristic lengths $s$ and $R$. The thickness of the sample is also small compared to $s$. However, the ratio $r/z$ is not necessarily negligible. In the cases when it is indeed negligible, the sample can be supposed planar instead of being volumetric (i.e., $r=0$) from a practical point of view; in this case, $t$ becomes a planar variable $(t_1,t_2)\in S$ (a square) and $m$ is actually a magnetization density.

The surface $Q_R$ of measurements is, most of the time, supposed to be a centered square of half-size $R$, but in some situations it might be convenient to consider only the data available on a centered disk $D_R$ of radius $R$. The microscope provides a map of the field on the whole surface: measurements are provided at many points of the surface and, from a practical point of view, one may assume that the field is known everywhere on $Q_R$. In addition to the presence of noise in the measurements, an important limitation is that, depending on the microscope technology, it is frequent that only one component of the field be measured.

The team has a long-standing collaboration with the Earth and Planetary Sciences Laboratory at MIT. They have a superconducting quantum interference device (SQUID). The sensor is a tiny vertical coil maintained at temperature close to 0 Kelvin, which provides it with superconducting characteristics. In order to maintain the sensor at very low temperature while the microscope operates in a room at normal temperature, the sensor is isolated behind a sapphire window. Though thin, this window enforces a measurement height $z$ such that $r\ll z\ll s$ and the sample can usually be supposed planar. Also, because the coil only allows to measure the field along its axis, the SQUID only provides measurements of $B_3$ and not the whole field.

Another type of microscopes consists in the quantum diamond microscopes (QDM). They use properties of special diamonds which, when properly excited with a laser and a microwave field, become luminescent in the presence of a magnetic field. From the difference of brightness of this luminescence under slightly different frequencies of the microwave field, one can recover the amplitude of a given component of the magnetic field. This mechanism is already in use to provide a magnetic microscope at Harvard University (Massachusetts, USA). We are collaborating with geoscientists of the Geophysics and Planetology Department of Cerege (CNRS, Aix-en-Provence) and physicists from ENS Paris Saclay to help them designing their own QDM. The promises of the QDM are manifolds, see also Section 7.3. First, they should allow measurements of the field at a height $z$ above the sample that is way smaller than what is permitted with the SQUID. This improves the spatial resolution and together with the conditioning of the inverse magnetization problem. However, this also imposes to model the sample as a 3-D object, as its thickness $r$ becomes usually comparable to $z$ in this case. Second, the technology of the QDM could make it possible, in principle, to measure the three components of the field instead of only one (see Section 7.4), which opens the way to new regularization techniques for the inverse problems. However, the measurement of a field with a QDM is more indirect than with a SQUID and this raises specific issues: in particular, it might be necessary to impose an external field on the sample, called a bias field, in order to resolve ambiguities when recovering the field from the brightness of the luminescence, and such a bias field can actually perturb the magnetic properties of the rock sample under study.

The issues raised by the inverse magnetization problem in the framework of magnetic microscopes such as SQUIDs or QDMs are numerous and we got several contributions on the subject over time. Particularly important for the full recovery of the magnetization are the silent sources, i.e., the magnetization that belong to the kernel of the direct operator, or in other terms, those magnetization that produce no field on the measurement area, see Section 3.2.2. We fully characterized such magnetizations in the thin-plate hypothesis (i.e., when $r$ is assumed to be 0), 40. Contrary to the full magnetization, the total net moment (i.e., the integral of the magnetization over the sample) is in principle a piece of information that could be retrieved from the measurements, since silent sources have a null net moment. In order to recover the total net moment, we proposed to use linear estimators and defined a bounded extremal problem to find good such linear estimators 3. However, additional hypotheses must be added in order to ensure the uniqueness of a solution for the full inversion problem. Such an hypothesis is provided when the support of the magnetization inside the sample is supposed to be sparse, e.g., composed of isolated points or 1-D curves. In this case, we proposed to use total variation regularization to solve the inverse problem 8, 50, see Section 3.2.3.

4.3 Inverse magnetization issues with the lunometer

Measurements of the remanent magnetic field of the Moon let geoscientists think that the Moon used to have a magnetic dynamo for some time, but the exact process that triggered and fed this dynamo is not yet understood, much less why it stopped. In particular, the Moon is too small to have a convecting dynamo like the Earth has. In order to address this question, our geoscientists colleagues at Cerege decided to systematically analyze the rock samples brought back from the Moon by Apollo missions.

The samples are kept inside bags with a protective atmosphere, and geophysicists are not allowed to open the bags, nor to take out the samples from NASA facilities. Moreover, the measurements must be performed with a passive device in order to ensure that the samples would not be altered by the measurements: in particular no cooling or heating is allowed, and neither is the use of anything producing a magnetic field like, e.g., motors. Finally, since the measurements must be performed directly at NASA, the instrument must be easy to take apart and to assemble once on site. The overall time devoted to measuring all samples is limited and each sample must be analyzed quickly (typically within a few minutes). For all these reasons, our colleagues from Cerege designed a specific magnetometer called the “lunometer”: this device provides measurements of the components of the magnetic field produced by the sample, at some discrete set of points located on disks belonging to three cylinders (see Figure 3). The goal was not to get a deep understanding of the magnetic properties of the studied samples with such a rudimentary instrument but rather to help selecting a few of them that seems really interesting to study in more details: this would be used to file a request to NASA to buy sub-samples of a few grams on which more instructive (though possibly destructive) experiments could be performed.

Show image description

The three orthogonal planes of the 3D canonical system of coordinates. Perpendicular to each of the planes, there is a cylinder and on each cylinder three sections are highlighted as circles drawn in black, blue and red. On each circle, regularly spaced dots illustrate the points where measurements are performed.

The collaboration with Cerege on this topic started in the framework of the MagLune ANR project whose overall objective was to devise models to explain how a dynamo phenomenon was possible on the Moon. Our contribution is to design methods to tell, from the measurements provided by the Lunometer, whether the remanent magnetization of the sample under study could be well modeled by a single magnetic dipole, and if so, what would be the position and magnetic moment of this dipole. To this end, we use ideas similar to those underlying the FindSources3D tool (see Sections 3.3, 3.5.4): we use rational approximation techniques to recover the position of the dipole; recovering the moment is then a rather simple linear problem. The rational approximation solver gives, for each circle of measurements, a partial information about the position of the dipole. These partial informations obtained on all nine circles must then be combined in order to recover the exact position. Theoretically speaking, the nine partial informations are redundant and the position could be obtained by several equivalent techniques. But in practice, due to the fact that the field is not truly generated by a single dipole, and also because of noise in the measurements and numerical errors in the rational approximation step, all methods do not show the same reliability when combining the partial results. We studied several approaches, testing them on synthetic examples, with more or less noise, in order to propose a good heuristic for the reconstruction of the position 73.

4.4 Shape identification of metallic objects

We started an academic collaboration with partners at LEAT (Laboratoire d'Électronique, Antennes, Télécommunications, Université Côte d'Azur – CNRS) on the topic of inverse scattering using frequency dependent measurements. As opposed to classical electromagnetic imaging where several spatially located sensors are used to identify the shape of an object by means of scattering data at a single frequency, a discrimination process between different metallic objects is here being sought for by means of a single, or a reduced number of sensors that operate on a whole frequency band. The spatial multiplicity and complexity of antenna sensors is here traded against a simpler architecture performing a frequency sweep.

The setting is shown on Figure 4. The total field $E_t$ is the sum of the incident field $E_{in}$ (here a plane wave) and scattered field $E_s$: at every point $X=(r,\theta ,\phi )$ in space we have $E_t=E_{in}+E_s$. A harmonic time dependency ($e^{i\omega t})$, is supposed for the incident wave, so that by linearity of Maxwell equations and after a transient state, the scattered field at the observation point $X_o$ is related to the emitted planar wave field at the emission point $X_e$ via the transfer function $H$: $E_s\left(X_o\right)=H(\omega ,X_o)E_{in}\left(X_e\right)$ (the dependency in $X_e$ is omitted in $H$ since the emission point is fixed). Under regularity conditions on the scatterer's boundary, the function $H$ can be shown to admit an analytic continuation into the complex left half-plane for the $s$ variable, away from a discrete set (with a possible accumulation point at infinity) where it admits poles. Thus, $H$ is a meromorphic function in the variable $s$. Its poles are called the resonating frequencies (resonances) of the scattering object. Recovering these resonating frequencies from frequency scattering measurement, that is measurements of $H$ at particular ${\omega }_j^{\text{'}}s$ (actually, $i{\omega }_j^{\text{'}}s$) is the primary objective of this project.

Show image description

A schematic view of the problem: an antenna (horizontal and on the left) emits a planar wave (propagating horizontally, towards the right and shown in red) with an emitted electric field $E$ perpendicular to it (represented by a vertical arrow on the figure). On the right, a sphere of radius $a$ reflects the wave. The reflected wave is radial to the sphere, directed towards its exterior, and is shown in blue. A reception antenna points towards the center of the sphere, and is located at a distance $r$ of the center of the sphere, while making an angle $\theta$ with respect to the horizontal. It measures the reflected electric field, which is perpendicular to it.

We started a study of the particular case when the scatterer is a spherical PEC (Perfectly Electric Conductor). In this case, Maxwell equations can be solved by means of expansions in series of vectorial spherical harmonics. We showed in particular that in this case $H$ admits a simple structure involving a meromorphic function with poles at zeroes of the spherical Hankel functions and their derivatives. Identification procedures, surprisingly close to the ones we developed in connection with amplifier stability analysis (Section 4.6), were studied to gain information about the resonating frequencies by means of a rational approximation of this function 91.

In order to perform the rational approximation (see Section 3.3), the behavior of $H$ outside the range of measured frequencies, specifically at high frequencies, has been studied for the particular case when the scatterer is a spherical PEC (Perfectly Electric Conductor). In this case, $H$ can be written as $H=H_O+H_C$, where $H_O$ and $H_C$ are respectively the optic and creeping wave parts. Their high-frequency behaviors are given by series expansions whose coefficients are identified when $X_o=X_e$. The asymptotics of $H_O$ is called the Luneberg-Kline expansion; its first terms were analytically computed (solving eikonal and transport equations), 31.

Numerical simulations showed that even though $H_C$ is negligible with respect to $H_O$ at high frequencies, it needs to be taken into account around the band of measured frequencies for the rational approximation. Furthermore, the physical interpretation of these two terms leads to consider that $H_C$ should carry more information about the scatterer and we want to investigate the conjecture that the poles of $H$ are those of $H_C$ hence that $H_O$ is analytic.

The rational approximation of the transfer function $H$ is performed with a least-squares substitute to multi-point Padé approximation: the approximant $R_{k,n}\left[H\right]=p_{k,n}/q_{k,n}$ is given by:

where $n$, $k$ and $N$ are integers such that $k+n+1\le N$ and ${\left(z_j\right)}_{j=1...N}$ is a collection of points at which $H$ is analytic. Here, ${𝒫}_n$ is the set of polynomials of degree less than $n$ and ${𝒫}_k^1$ is the set of monic polynomials of degree $k$. An analog of the Nuttal-Pommerenke theorem for a least square version of classical Padé approximants was obtained in 31, where the values $(p-Hq)\left(z_j\right)$ in (1) get replaced by the $j$-th coefficients of the Taylor expansion of $(p-Hq)$ at a given point in the domain of analyticity of $H$. It says that if $H$ is holomorphic and single-valued on $ℂ$, except perhaps on a polar set, then $p/q$ converges to $H$ in capacity as $k$, $n$ go to infinity as well as $N$, in such a way that $n/k$ remains bounded and $N\le C(k+n)$ for some constant $C>0$; convergence in capacity means that for each $\epsilon >0$, the capacity of the set where the pointwise error is bigger than $\epsilon$ goes to zero.

Dwelling on this work, we recently showed under similar conditions on $k$, $n$, $N$, and provided that the interpolation points $z_i$ remain in a bounded set, that convergence in capacity still holds for the solutions to (1); this is a least square analog of Wallin's theorem in multi-point Padé approximation.

This result is interesting as it entails that poles of $H$ can indeed be retrieved as limit points of certain poles of $p_{k,n}/q_{k,n}$, while explaining the chaotic behavior of other, so-called spurious poles that wander about the domain of analyticity. We plan to investigate other PEC scatterers.

4.5 Inverse problems in orthopedic surgery

Apart from more classical medical imaging domains, inverse problems find a rather surprising application in the field of orthopedics, see 60, 70, 74 and Section 8.2.

We are concerned, in particular, with a hip prosthetic surgery when an insertion of an acetabular cup (AC) implant into a bone by press-fit is performed with the use of an instrumented hammer. Such a hammer is equipped with a sensor capable to measure impact momentum (force) and hence yield important information about the bone quality and the stability of an AC implant. These are, indeed, crucial pieces of information to have in real-time during a surgery. On the one hand, if the achieved bone-implant contact area is not sufficiently large, osseointegration may fail eventually leading to an aseptic loosening of the implant and a necessity of a revision surgery. On the other hand, if the insertion of the implant is too deep, the generated stresses may induce fractures or bone tissue necrosis.

The mathematical side of the process is far from trivial. First of all, contact mechanics is a highly nonlinear problem due to geometrical constraint on the solution. Already a basic problem of an elastic body on a rigid foundation is a free-boundary problem with an unknown effective contact surface which is characterized by the so-called Signorini conditions, nonlinear constraints of Karush-Kuhn-Tucker type involving stress and displacements. In the present case, several coupled problems have to be solved since regions corresponding to the bone, the implant and the hammer all possess different material properties. Moreover, the deformations cannot be considered small, consequently, a hypo-elastic description is more appropriate than that of linear elasticity. In such a formulation, a rate relation between Cauchy stress and strain tensors replaces a linear stress-strain constitutive law, hence its integration induces additional non-linearity. Finally, the bone is a porous multi-scale medium, and appropriate homogenization model should be deduced, with adequate parameter fitting.

For tackling inverse problems (e.g., that of determining material parameters), the direct formulation has to be solved in such an effective way that iterative approaches are not prohibitively expensive. This motivates exploration of model-order reduction strategies that would, in particular, allow efficient integration of the system in time.

4.6 Stability and design of active devices

Through contacts with CNES (Toulouse) and UPV (Bilbao), the team got involved in the design of amplifiers which, unlike filters, are active devices. A prominent issue here is stability. Twenty years ago, it was not possible to simulate unstable responses, and only after building a device could one detect instability. The advent of so-called harmonic balance techniques, which compute steady state responses of linear elements in the frequency domain and look for a periodic state in the time domain of a network connecting these linear elements via static non-linearities made it possible to compute the harmonic response of a (possibly nonlinear and unstable) device 87. This has had tremendous impact on design, and there is a growing demand for software analyzers. In this connection, there are two types of stability involved. The first is stability of a fixed point around which the linearized transfer function accounts for small signal amplification. The second is stability of a limit cycle which is reached when the input signal is no longer small and truly nonlinear amplification is attained (e.g., because of saturation).

Initial applications by the team have been concerned with the first type of stability, and emphasis was put on defining and extracting the “unstable part” of the response. We showed that under realistic dissipativity assumptions at high frequency for the building blocks of the circuit, the linearized transfer functions are meromorphic in the complex frequency variable $s$, with at most finitely many unstable poles in the right half-plane 2. Dwelling on the unstable/stable decomposition in Hardy Spaces, we developed a procedure to assess the stability or instability of the transfer functions at hand, from their evaluation on a finite frequency grid 10, that was further improved in 55 to address the design of oscillators. This has resulted in the development of a software library called Pisa (see Section 3.5.1), aiming at making these techniques available to practitioners. Extending this methodology to the strong signal case, where linearization is considered around a periodic trajectory, is considerably more difficult and has received much attention by the team in recent years. The exponential stability of the high frequency limit of a circuit was established in 5, implying that there are at most finitely many unstable poles and no other unstable singularity for the monodromy operator around the cycle. Furthermore, the links between the monodromy operator and the (operator-valued) “harmonic transfer function” (HTF) of the linearized system along the trajectory were brought to light in 24: the system is exponentially stable if and only if its HTF is bounded an analytic in a right half-plane of $ℂ$ of the form $\{\Re z>-\epsilon \}$ for some $\epsilon >0$. We deal here with input-output system of the form:

where ${\tau }_1<\cdots <{\tau }_N$ are positive delays and $D_1\left(t\right),...,D_N\left(t\right)$ real $d\times d$ matrices depending periodically on time $t$, while $u$ is the ${ℝ}^d$-valued input and $y\left(t\right)$ the ${ℝ}^d$-valued output (complex coefficients can of course be handled in the same way). The HTF, that generalizes the usual transfer function of linear constant systems, is an analytic function of the Laplace-Fourier variable which is valued in the space of operators on $L^2([0,T),{ℂ}^d)$ ($T$ is the period of the system). It can be defined as follows: if a periodic linear control system at rest is fed from initial instant $t_i=-\infty$ with an input signal $v\left(t\right)e^{xt}{e}^{i\nu t}$ where $v$ is $T$-periodic and $x>0$ is large enough, then the output is of the form $𝐇(x+i\nu )\left(v\right)e^{xt}{e}^{i\nu t}$, where $𝐇(x+i\nu )$ is the harmonic transfer evaluated at $x+i\nu$ (an operator) and $𝐇(x+i\nu )\left(v\right)$ is the $T$-periodic function which is the image of the $T$-periodic function $v$ under $𝐇(x+i\nu )$. That is to say, an exponentially modulated input wave of frequency $\nu$ carried by a $T$-periodic signal is mapped to an output of the same form, and the HTF maps the input carrier to the output carrier. For stable systems, this is asymptotically true if the initial time is a finite instant; otherwise, the unstable transient will prevent this mathematical solution from being physical. Other definitions are given in 24.

We were able to construct a simple nonlinear circuit whose linearization around a periodic trajectory has a spectrum containing a whole circle; we currently investigate whether the singularities of the Fourier coefficients of the HTF (that are themselves analytic functions) also contain that circle. Indeed, just like a series of functions may diverge even though the summands are smooth, it is a priori possible that the HTF has a singularity at a point whereas its Fourier coefficients do not. Note that these coefficients are all one can estimate by harmonic balance techniques, and therefore the above question is of great practical relevance. We also investigate the relation between our stability criterion (that the HTF should be bounded and analytic on a “vertical” half-plane containing the origin), and the weaker requirement that the HTF exists pointwise in a half-plane. Connections with the representation of Volterra equations with jumps in the kernel are also a motivation for such a study, see 16.

4.7 Tools for numerically guaranteed computations

The overall and long-term goal is to enhance the quality of numerical computations. This includes developing algorithms whose convergence is proved not only when assuming that the numerical computations are performed in exact real or complex arithmetic, but rather when really accounting for the fact that the computations are performed with an inexact arithmetic (usually floating-point arithmetic). A numerical result alone is of little interest if no rigorous bound is provided together with it, in order to ensure that the real theoretical result is proved to be not to far from the computed result.

A specific way of contributing to this objective is to develop efficient numerical implementations of mathematical functions with rigorous bounds. We do sometimes provide such implementations. The software tool Sollya (see Section 3.5.6), developed together with C. Lauter (University of Texas at El Paso, UTEP) is also an achievement of the team in that respect. This tool intends to provide an interactive environment for performing numerically rigorous computations. Sollya comes as a standalone tool and also as a C library that allows one to benefit from all the features of the tool in C programs.

4.8 Imaging and modeling ancient materials

This is an additional activity of the team, linked to image classification in archaeology, pursued in collaboration with the partners listed in Section 8.3.

The pottery style is classically used as the main cultural marker within Neolithic studies. Archaeological analyses focus on pottery technology, and particularly on the first stages of pottery manufacturing processes. These stages are the most demonstrative for identifying the technical traditions, as they are considered as crucial in apprenticeship processes. Until now, the identification of pottery manufacturing methods was based on macro-traces analysis, i.e., surface topography, breaks and discontinuities indicating the type of elements (coils, slabs, ...) and the way they were put together for building the pots. Overcoming the limitations inherent to the macroscopic pottery examination requires a complete access to the internal structure of the pots. Micro-computed tomography can be used for exploring ancient materials micro-structure. This non-invasive method provides quantitative data for a big set of proxies and is perfectly adapted to the analysis of cultural heritage materials.

The main challenge of our current analyses aims to overcome the lack of existing protocols to apply in order to quantify observations. In order to characterize the manufacturing sequences, the mapping of the paste variability (distribution and composition of temper) and the discontinuities linked to different classes of pores, fabrics and/or organic inclusions appears promising. The totality of the acquired images composes a set of 2-D and 3-D data at different resolutions and with specific physical characteristics related to each acquisition modality. Specific shape recognition methods need to be developed by application of robust imaging techniques and 3-D-shapes recognition algorithms.

We devised in 54 a method to isolate pores from the 3-D data volumes in binary 3-D images, to which we apply a process named Hough transform (derived from Radon transform). This method, of which the generalization from 2-D to 3-D is quite recent, allows us to evaluate the presence of parallel lines going through the pores. The quantity of such lines and their parallelism furnish good indicators of the “coiling” manufacturing, that they allow to distinguish from the other “spiral patchwork” technique, in particular. Other possibilities of investigation will be analyzed, such as deep learning techniques.

6.1 Awards

Dmitry Ponomarev was awarded a runner-up poster prize at the 4th IMA Conference on Inverse Problems in Bath (September), see Section 9.1.4.

6.2 Contract of Objectives, Means and Performance

Note : Readers are advised that the Institute does not endorse the text in the “Highlights of the year” section, which is the sole responsibility of the team leader.

Inria management has imposed on the institute a new Contract of Objectives, Means and Performance (COMP) with the French government, for the period 2024–2028. It presages major changes for Inria, regarding both its missions and the way it operates. These changes, whose precise nature and impact on the staff are unclear, should become effective as soon as 2025 but have not been the subject of any consultation, and inasmuch as the collaboration of Inria's staff is necessary to turn this disruption into a successful change, we are concerned that the top management has remained deaf to several votes and petitions opposing these policies.

The multiplication of new missions and priorities, particularly those related to the “program agency” or oriented towards defense applications, pushes the research carried out at Inria in the background. The constraints induced by this COMP will restrain the independence of scientists and teams, as well as their freedom to select research topics and collaborators.

Our main concern with the COMP lies with the following items:

• Placement of Inria in a “zone à régime restrictif” (ZRR).
• Restriction of international and industrial collaborations to partners chosen by the institute's management, with no clear indication of the rules.
• Individual financial incentives for researchers involved in strategic partnerships, whose topics are steered by the program agency.
• Priority given to “dual” research with both military and civilian applications, materialized by tighter links with the Ministry of Defense.

7.1 Field extrapolation for inverse magnetization problem

Participants: Dmytro Dmytrenko, Juliette Leblond, Mubasharah Khalid Omer, Dmitry Ponomarev.

A simple fact that more measured data must yield better results manifests itself as a slower decay rate of singular values of the forward operator $A$ (introduced in Section 3.1.3) for a larger measurement region. The same phenomenon happens when the abundance of data is not only quantitative but qualitative. Namely, for the inverse magnetization problem (see Section 4.2), reconstruction results demonstrate obvious improvements when measurements of all three components of the magnetic field are used instead of only one. Because of the harmonicity of field, the knowledge of only one (normal) component of the field actually determines two other (tangential) components. However, this operation (Riesz transformation) is stable when the normal component of the field is known on the entire horizontal plane. Since in practice the actual measurements of the field data are limited to a portion of the plane, we are thus led to consider an extrapolation problem. This problem is well-defined in the ideal case of knowledge of exact data, i.e., one can uniquely construct the field defined on the whole plane which agrees with the given data on its portion. In case of perturbed data, the problem lacks stability but can nevertheless be effectively addressed. To this effect, we have pursued the development of two approaches.

The first approach is fairly straightforward and is based on the regularized deconvolution of the field in order to find an auxiliary magnetization-dependent function on ${ℝ}^2$ (see Section 3.1.2). Re-application of the Poisson transform to this function allows us to generate the field on the whole plane at height $x_3=h$ thus obtaining an extrapolant. This approach can be efficiently implemented when an appropriate finite-dimensional subspace is chosen. A particular choice that was investigated was the span of spherical harmonics mapped to the horizontal plane by a suitable Kelvin transform. Since Kelvin transform preserves harmonicity and orthogonality, these functions are naturally restrictions of harmonic functions and application of the Poisson operator to them is computationally cheap; harmonic extension of spherical harmonics can be obtained without any convolution, these are simply solid harmonics. Other choices of basis functions can be efficiently arranged products of orthogonal systems on $L^2\left(ℝ\right)$ such as Laguerre or Malmquist-Takenaka functions. This extrapolation approach is developed for planar $L^2$ magnetization distributions (see Section 3.1.2), but it does not profit from a usually known compact support of the magnetization. On the other hand, since a compact support is not a requirement, it can also be an advantage allowing to use this approach for three-dimensional magnetizations. Indeed, a volumetric sample is easily shown to be equivalent to a planar magnetization distribution (supported on the entire plane) at the height corresponding to the distance between the measurement plane and the original sample. This is the topic of the PhD thesis of Mubasharah Khalid Omer .

The second extrapolation approach 18, 19, 23, 27 assumes planar magnetization with square-integrable divergence of its tangential components and relies on the knowledge of its support (which is the case in practice). Namely, it is assumed that the planar measurement region covers the magnetization support in the underlying plane. This method relies on the representation of the field as sum of two convolution operators acting on a combination of tangential components of magnetization (actually, the divergence thereof) and the normal one, respectively. This involves kernel functions essentially coinciding with Poisson kernel and its vertical derivative. The extrapolated field is generated from the same relation once the two mentioned magnetization-dependent quantities are found. To this effect, we have developed a so-called double-spectral approach (in the spirit of Section 3.1.3). The approach consists in 2 stages. First, an additional (data-independent) problem is solved to generate an auxiliary kernel function whose convolution swaps (to a certain degree) the action of two other convolution operators. This can be achieved by the spectral decomposition of a truncated Poisson operator, i.e., explicit computation of its eigenvalues and eigenfunctions. At the second stage, the problem is immersed in a higher dimension and the self-adjoint system of integral equations is solved again by a spectral expansion. The use of spectral decomposition for the solution of the auxiliary problem at the first stage allows to efficiently characterize the null-space of the matrix operator occurring at the second stage. The use of another spectral method at the second stage (i.e., using eigen-elements of the matrix operator as a solution basis) is convenient for de-noising possibly perturbed field measurements by means of truncation of higher-order modes. The number of basis elements at both stages thus serve as regularization parameters which have to be adjusted according to the level of noise in measurements or to the accuracy of numerical computations of auxiliary quantities. With this respect, it should be mentioned that since the field extrapolation problem is subject to more constraints than that for reconstruction of the full magnetization, the additional relations (such as vanishing of total integral of $B_3$ over ${ℝ}^2$) can be employed for choosing regularization parameters.

7.2 Net moment estimation

Participants: Sylvain Chevillard, Juliette Leblond, Jean-Paul Marmorat, Anass Yousfi.

We continued the work started in the past years to establish formulas for the integrals of the form

where $P$ is a polynomial and the domain of integration is either the square $Q_R$ (and actually more generally, any rectangle) or the disk $D_R$ of radius $R$ (the notations are those of Section 4.2, see especially Figure 2). This is the topic of the PhD thesis of Anass Yousfi .

On the one hand, in the case of the disk, exact formulas are not available but we had obtained last year asymptotic expansions with respect to variable $R$ as it grows large. Due to unexpected personal difficulties, the report that we were working on has been delayed and is still on-going work. However, progresses were achieved by simplifying some proofs, while formulas have been implemented and successfully compute the sought asymptotic expansion.

On the other hand, in the case of the square, we obtained exact formulas, not only for $B_3$ but also for the other components of the field. More specifically, denoting by abuse of notations $B_j$ the (forward) operator from $L^2\left(S\right)$ to $L^2\left(Q_R\right)$ that maps a magnetization to the corresponding component of the field ($j\in \{1,2,3\}$), its adjoint operator $B_j^{☆}$ plays an important role for the estimation of the net moment from measured data. Recursion formulas allow one to compute an explicit form for $B_j^{☆}\left(P\right)$ for any polynomial $P$ 22. Our goal is now to see how this can be used for the evaluation of the solution to the extremal problem presented in 3, by expanding it on the basis of the eigenfunctions of the Laplace operator as suggested in the above cited article. We did some preliminary experiments in that direction by comparing the efficiency of two methods for the evaluation of integrals of the form ${\int \int }_{Q_R}sin\left(k_1{x}_1\right)sin\left(k_2{x}_2\right)B_3\left(x\right)\text{d}x$ (when $k_1$ and $k_2$ go large): either applying classical quadrature rules or approximate the factor $sin\left(k_1{x}_1\right)sin\left(k_2{x}_2\right)$ by a polynomial $P(x_1,x_2)$ and apply the formulas described above for the operator $B_3^{☆}$. Another kind of quadrature formulas, specifically tuned for oscillatory integrals also retained our attention: these are Filon quadrature rules 72. It is not clear yet which combination of these three strategies is optimal depending on the value of $k_1$ and $k_2$. Progresses on that respect are expected.

7.3 QDM magnetometry

Participants: Sylvain Chevillard, Juliette Leblond, Dmitry Ponomarev, Anass Yousfi.

Compared to SQUID microscopy (which provided the team with a concrete context for developing methodology for the inverse magnetization problem), Quantum-Diamond Magnetometric (QDM) devices are capable of producing high-resolution field maps, potentially of more than one component, and extremely close to a magnetic sample 62, 82, 88. One of such devices is currently assembled in Cerege that the team already has a collaboration with (see Sections 4.2, 4.3).

The measurement principle of QDM can be roughly described as follows. A point defect in the crystalline lattice of diamond is introduced and consists of a nitrogen-vacancy (NV) pair. Quantum structure of this NV center leads to its crucial magneto-optical properties. The NV center is a two-spin system (2 out of its 6 electrons with a spin $1/2$ are unpaired) whose triplet state (i.e., that with the total spin equal to 1) is of particular interest. Denoting ${𝚖}_{𝚜}$ the projection of the total spin of the NV center onto its axis, the following is observed. When excited by laser, NV states with $\left|{𝚖}_s\right|=1$ and ${𝚖}_s=0$ relax back to the ground state differently causing the drop in signal when diamond's illumination intensity is measured for different frequencies. Furthermore, when an exterior magnetic field is added, it lifts degeneracy of the states ${𝚖}_s=\pm 1$ according to the Zeeman effect and thus causes the splitting of this intensity peak (drop) in two. Difference of frequencies of these two peaks is proportional to the magnitude of projection of the field onto NV access which gives a possibility of a magnetic measurement. Since a NV center can be, in general, oriented along 4 different directions (according to the diamond crystalline lattice symmetry), instead of 2, 8 peaks are observed experimentally, and reconstruction of the magnetic field through absolute value of its projections on those axes is not a trivial task.

Mathematical contribution here can lead to a better design of experimental set-ups. As a first step, we have taken up an issue of improving a quantitative estimate of the Zeeman effect relating the distance between peaks to the field projection magnitude beyond the linear order and estimating the error. This is important because the transverse (with respect to the NV axis) component of the field may not necessarily be small.

7.4 Inverse magnetization problem with 3-component field measurements

Participants: Laurent Baratchart, Sylvain Chevillard, Juliette Leblond, Dmitry Ponomarev.

Since missing components of the magnetic field can be either obtained by extrapolation (see Section 7.1), or by QDM measurements (see Section 7.3), it is of major interest to extend the methodology of previously developed approaches (dealing with different aspects of inverse magnetization problems) to take into account availability of this new type of information. This concerns, for example, reconstruction of net moments (see Sections 3.2.2, 4.2, 7.2) by means of asymptotic estimates 36, 42, 81 and through solution of an auxiliary bounded extremal problem 3. A preliminary work in this direction has already started. As already mentioned in Section 7.1, the results for the magnetization reconstruction (its minimal $L^2$-norm representative, among equivalent magnetizations producing the same field data) with $L^2$ constraints are shown to be significantly better when all three components of the field are used. Numerical implementation, however, depends on how well the magnetization is represented in a chosen finite-dimensional basis. In particular, we have used suitably ordered tensor products of one-dimensional modified Fourier basis functions (eigenfunction of Laplace operator with zero Neumann boundary data) which are great for smooth magnetizations, but not so much for “grainy” ones. For the universal choice, one perhaps should proceed with a wavelet-type basis which can represent multi-scale features and consider its adaptation taking into account the exact form of “silent” magnetic sources.

7.5 Rational and meromorphic approximation

Participants: Laurent Baratchart.

In a joint research effort with M. Yattselev from Indiana University at Indianapolis, we completed this year a long term project that started years ago as a collaboration with the late H. Stahl. It deals with the behavior of rational approximants of given degree $N$ to a holomorphic function $f$, on a compact subset $O$ of the domain of analyticity bounded by Jordan curve $\gamma$, under the assumption that $f$ can be continued analytically to the exterior of $\gamma$ (possibly in a multiple-valued manner) except perhaps over a set of logarithmic capacity zero of the extended plane.

If we let $R_N$ indicate rational functions of degree $N$ with no pole in $O$, a best rational approximant is a member $r_N$ of $R_N$ minimizing $\parallel f-r_N{\parallel }_{L^{\infty }\left(\gamma \right)}$; we say that a sequence $r_N$ of rational approximants, with $r_N\in R_N$, is asymptotically optimal in the $N$-th root sense if ${lim\; sup}_N{\parallel f-r_N\parallel }_{L^{\infty }\left(\gamma \right)}^{1/N}$ is as small as possible (see also Section 3.3).

The behavior of rational approximants is companion to the one of meromorphic approximants to $f$ on $\gamma$. To define the latter, let $R_N$ indicate rational functions of degree $N$ with no pole in $O$ and ${𝒜}_{\gamma }$ be the space of bounded holomorphic functions in the exterior of $\gamma$, that extend continuously up to $\gamma$. For our present purpose, a meromorphic function is the trace on $\gamma$ of a member of $M_N:={𝒜}_{\gamma }+R_N$, and a best meromorphic approximant solves the following extremal problem:

($E_N$)  Let $N\ge 0$ be an integer, and $f$ as above; to find a function $g_N\in M_N$ such that $g_N-f$ is of minimal norm in $L^{\infty }\left(\gamma \right)$.

That $\left(E_N\right)$ is well-posed depends on the fact that $f$ extends analytically across $\gamma$ by assumption. When $\gamma$ is rectifiable, the unique solution is given by AAK theory (named after Adamjan, Arov and Krein) in $\overline{ℂ}\setminus O$, which connects the spectral decomposition of Hankel operators with best approximation 78. Again, we say that a sequence $g_N$ of meromorphic approximants, with $g_N\in R_N$, is asymptotically optimal in the $N$-th root sense if ${lim\; sup}_N{\parallel f-g_N\parallel }_{L^{\infty }\left(\gamma \right)}^{1/N}$ is as small as possible. It follows from the work in 63, 89 that, for $f$ as above, ${lim\; sup}_N{\parallel f-r_N\parallel }_{L^{\infty }\left(\gamma \right)}^{1/N}$ is a true limit whenever $r_N$ is asymptotically optimal in the $N$-th root sense, and this limit is equal to $exp\{-2/\mathrm{cap}(O,{𝒦}_f)\}$ where ${𝒦}_f$ is the compact set minimizing the capacity of the condenser $(O,{𝒦}_f)$ outside of which the extension of $f$ is single-valued. Moreover, the same holds for meromorphic approximants that are asymptotically optimal in the $N$-th root sense, by 76. The existence and uniqueness of ${𝒦}_f$, up to a polar set, is a consequence of 84. If the extension of $f$ has no branch-point then ${𝒦}_f=\varnothing$ and the convergence of rational or meromorphic approximants that are asymptotically optimal in the $N$-th root sense is faster than geometric, otherwise $\mathrm{cap}(O,{𝒦}_f)>0$ and the convergence is geometric.

Now, if the extension of $f$ outside $\gamma$ has finitely many branches; that is, if it is analytic and single valued on a compact Riemann surface $𝒮$ (except possibly on a polar subset of $𝒮$), then we were able to show that rational or meromorphic approximants that are asymptotically optimal in the $N$-th root sense converge in capacity in the complement of ${𝒦}_f$, meaning that for any $\epsilon >0$ the set of points where the error is larger than $\epsilon$ has a logarithmic capacity that tends to zero when $N$ goes to infinity. Moreover, if $𝒮$ is a sphere (equivalently: if the continuation of $f$ has no branch-point), then this convergence is faster than geometric and there is a sequence of approximants, converging at least as fast, whose poles coalesce to the polar singular set of $f$. And if $𝒮$ has positive genus (equivalently: if the continuation of $f$ has branch-points), then the normalized counting measures of the poles of the approximants converge weak-star, as $N$ tends to infinity, to the equilibrium distribution on the plate ${𝒦}_f$ of the condenser $(O,{𝒦}_f)$; see 25. Interpolation techniques and complex orthogonal polynomials, that form the technical core of the work in 63 are of no help in our case, as it is not even known if interpolation takes place for best approximants in the sup norm. The proof proceeds from completely different ideas, reducing to AAK-type of approximants and analyzing sub-level sets of the logarithmic error function on the underlying Riemann surface; we refer the reader to 25 for details.

Let us also mention that convergence in capacity of least square substitutes to multi-point Padé approximants (1) to functions with polar singular set on $ℂ$ (see Section 4.4 for a precise description) was established by elaborating on the techniques developed in 31 and reworking normalization issues. An article is being written to report on this result.

7.6 Schur functions minimization under Nevanlinna-Pick constraints

Participants: Martine Olivi, Fabien Seyfert.

In the past, Factas has devoted a great deal of effort to the problem of impedance matching in filter synthesis. Filter synthesis is usually performed under the hypothesis that both ports of the filter are loaded on a constant, resistive load (usually 50 $\Omega$). In complex systems, however, filters are cascaded with other devices and end up being loaded on a non purely resistive, frequency varying, load. The load thus varies with frequency by construction, and is not purely resistive either. Likewise, in an emitter-receiver, the antenna is followed by a filter; however, the antenna is far from being purely resistive on the whole passband, and a mismatch between the antenna and the filter then causes irremediable power losses, both for emission and transmission. Our goal has therefore been to develop methods for filter synthesis that allow us to match varying loads on specific frequency bands, while enforcing rejection properties away from the passband. The problem is highly nontrivial because matching requires the response of the filter, which is holomorphic in the right half-plane (since the filter is a stable device), to approximate a function which is holomorphic in the left half-plane.

An approach for solving this issue amounts to consider the following optimization problem: find a Schur function with minimum supremum norm in a frequency band, which satisfies a set of Nevanlinna-Pick interpolation constraints. A convex approach to this problem was developed by D. Martínez Martínez and G. Bose during their PhD theses. The preprint 26 synthesizes their work and generalizes it to the case of a mixed Nevanlinna-Pick problem, with interpolation points both inside and on the boundary of the analyticity domain.

8.1 International research visitors

8.2 National initiatives

8.3 Regional initiatives

9.1 Promoting scientific activities

9.2 Teaching - Supervision - Juries

9.3 Popularization

10.1 Major publications

• 1 articleB.Bilal Atfeh, L.Laurent Baratchart, J.Juliette Leblond and J. R.Jonathan R. Partington. Bounded extremal and Cauchy-Laplace problems on the sphere and shell.J. Fourier Anal. Appl.162Published online Nov. 20092010, 177--203URL: http://dx.doi.org/10.1007/s00041-009-9110-0back to text
• 2 articleL.Laurent Baratchart, S.Sylvain Chevillard, A.Adam Cooman, M.Martine Olivi and F.Fabien Seyfert. Linearized Active Circuits: Transfer Functions and Stability.Mathematics in Engineering45October 2021, 1-18HALDOIback to text
• 3 articleL.Laurent Baratchart, S.Sylvain Chevillard, D. P.Douglas P. Hardin, J.Juliette Leblond, E. A.Eduardo Andrade Lima and J.-P.Jean-Paul Marmorat. Magnetic moment estimation and bounded extremal problems.Inverse Problems and Imaging 131February 2019, 29HALDOIback to textback to textback to text
• 4 articleL.Laurent Baratchart, S.Sylvain Chevillard and T.Tao Qian. Minimax principle and lower bounds in $H^2$-rational approximation.Journal of Approximation Theory2062015, 17--47back to text
• 5 articleL.Laurent Baratchart, S.Sébastien Fueyo, G.Gilles Lebeau and J.-B.Jean-Baptiste Pomet. Sufficient Stability Conditions for Time-varying Networks of Telegrapher's Equations or Difference Delay Equations.SIAM Journal on Mathematical Analysis532021, 1831–1856HALDOIback to text
• 6 articleL.Laurent Baratchart, J.Juliette Leblond, S.Stéphane Rigat and E.Emmanuel Russ. Hardy spaces of the conjugate Beltrami equation.Journal of Functional Analysis25922010, 384-427URL: http://dx.doi.org/10.1016/j.jfa.2010.04.004back to textback to text
• 7 articleL.Laurent Baratchart, H.Herbert Stahl and M.Maxim Yattselev. Weighted Extremal Domains and Best Rational Approximation.Advances in Mathematics2292012, 357-407URL: http://hal.inria.fr/hal-00665834back to text
• 8 articleL.Laurent Baratchart, C.Cristobal Villalobos Guillén, D.Doug Hardin, M.Michael Northington and E.Edward Saff. Inverse Potential Problems for Divergence of Measures with Total Variation Regularization.Foundations of Computational Mathematics202020, 1273-1307HALDOIback to textback to text
• 9 articleM.Maureen Clerc, J.Juliette Leblond, J.-P.Jean-Paul Marmorat and T.Théodore Papadopoulo. Source localization using rational approximation on plane sections.Inverse Problems285May 2012, Article number 055018 - 24 pagesHALDOIback to textback to text
• 10 articleA.Adam Cooman, F.Fabien Seyfert, M.Martine Olivi, S.Sylvain Chevillard and L.Laurent Baratchart. Model-Free Closed-Loop Stability Analysis: A Linear Functional Approach.IEEE Transactions on Microwave Theory and Techniqueshttps://arxiv.org/abs/1610.032352017HALback to text
• 11 inbookJ.Juliette Leblond and D.Dmitry Ponomarev. On some extremal problems for analytic functions with constraints on real or imaginary parts.Advances in Complex Analysis and Operator Theory: Festschrift in Honor of Daniel Alpay's 60th BirthdayBirkhauser2017, 219-236HALDOIback to text
• 12 articleM.Martine Olivi, F.Fabien Seyfert and J.-P.Jean-Paul Marmorat. Identification of microwave filters by analytic and rational H2 approximation.Automatica492January 2013, 317-325HALDOIback to text

10.2 Publications of the year