On the travel time tomography problem in 3D

TL;DR

This paper addresses the 3D travel time tomography problem with limited data, introducing a stable numerical method based on Fourier series and Carleman estimates to ensure global convergence.

Contribution

It develops a novel globally convergent numerical approach for 3D travel time tomography using Fourier series and Carleman estimates, handling non-overdetermined data.

Findings

Lipschitz stability estimate established

A globally convergent numerical method constructed

Effective handling of non-overdetermined data

Abstract

Numerical issues for the 3D travel time tomography problem with non-overdetemined data are considered. Truncated Fourier series with respect to a special orthonormal basis of functions depending on the source position is used. In addition, truncated trigonometric Fourier series with respect to two out of three spatial variables is used. First, the Lipschitz stability estimate is obtained. Next, a globally convergent numerical method is constructed using a Carleman estimate for an integral operator.

Equations

$$G=\left\{\mathbf{u}=\left(x,y,z\right):0<x,y<\pi ,z\in \left(B_1,B_2\right)\right\}.$$

$${\Gamma }_l=\left\{0<x,y<\pi ,z=B_1\right\},{\Gamma }_u=\left\{0<x,y<\pi ,z=B_2\right\},$$

$${\partial }_{1}G=\partial G╲\left({\Gamma }_l\cup {\Gamma }_u\right).$$

$$c\in C^2\left(\overline{G}\right),c\in C\left({\mathbb{R}}^3\right),$$

$$c\left(\mathbf{u}\right)\ge c_0,\forall \mathbf{u}\in {\mathbb{R}}^3,$$

$${\partial }_{z}c\left(\mathbf{u}\right)\ge 0,\mathbf{u}\in G,$$

$$c\left(\mathbf{u}\right)=1\text{ for }z<B_1.$$

$$dt=\sqrt{c\left(\mathbf{u}\right)}\left|d\mathbf{u}\right|,|d\mathbf{u}|=\sqrt{(dx{)}^2+(dy{)}^2+(dz{)}^2}.$$

$$t\left(\mathbf{u},{\mathbf{u}}_0\right)=\underset{\Phi \left(\mathbf{u},{\mathbf{u}}_0\right)}{\int }\sqrt{c\left(\xi \right)}d\sigma .$$

$$t_z^2+t_x^2+t_y^2=c\left(\mathbf{u}\right),$$

$$t\left(\mathbf{u},{\mathbf{u}}_0\right)=O\left(\left|\mathbf{u}-{\mathbf{u}}_0\right|\right),\mathbf{u}\to {\mathbf{u}}_0.$$

$${\mathbf{u}}_a\in I=\left\{\left(x,y,z\right):x=a\in \left[0,\pi \right],y=\pi /2,z=0\right\}.$$

$$t\left(\mathbf{u},a\right)=p\left(\mathbf{u},a\right),\forall \mathbf{u}\in \partial G╲{\Gamma }_l,\forall a\in \left(0,\pi \right).$$

$$t\left(\mathbf{u},a\right)=\sqrt{{\left(x-a\right)}^2+{\left(y-\pi /2\right)}^2+B_1^2}\text{ for }\mathbf{u}=\left(x,y,B_1\right)\in {\Gamma }_l,a\in \left[0,\pi \right].$$

$${\left\{{\phi }_n\left(a\right)\right\}}_{n=0}^{\infty }={\left\{P_n\left(a\right)e^a\right\}}_{n=0}^{\infty }$$

$$D_{N}f=\sum _{m=0}^{N-1}{f}_m{\phi }_m\left(a\right),f_m=\left(f,{\phi }_m\right).$$

$$\frac{d}{da}\left(D_{N}f\right)=\sum _{n=0}^{N-1}{f}_m{\phi }_m^{\prime }\left(a\right).$$

$$\left(D_{N}f\right)\in C^{\infty }\left(0,a\right),\frac{d}{da}\left(D_{N}f\right)\in L_{2,N}\left(0,\pi \right),\forall f\in L_2\left(0,a\right).$$

$$y_n^N=\left(\frac{d}{da}\left(D_{N}f\right),{\phi }_n\right),y^N={\left(y_0^N,...,y_{N-1}^N\right)}^T,z_N={\left(f_0,...,f_{N-1}\right)}^T.$$

$$z_N=A_N^{-1}\left(y^N\right).$$

$$\frac{d}{da}\left(D_{N}f\right)\ne D_N\left(f^{\prime }\right).$$

$$\Phi \left(\mathbf{u},a\right)=\left\{\left(x,y,z\right):x=x\left(z,a\right),y=y\left(z,a\right),z\in \left(0,Z\left(\mathbf{u},a\right)\right)\right\},$$

$$t_z^2\left(\mathbf{u},a\right)=\underset{B_1}{\overset{z}{\int }}{c}_z\left(x\left(s,a\right),y\left(s,a\right),s\right)ds+t_z^2\left({\mathbf{v}}_a,a\right).$$

$$t_z^2\left({\mathbf{v}}_a,a\right)=\frac{B_1^2}{{\left(x_a-a\right)}^2+{\left(y-\pi /2\right)}^2+B_1^2}\ge \frac{B_1^2}{B_1^2+5{\pi }^2/4}.$$

$$t_z\left(\mathbf{u},a\right)\ge \frac{B_1}{\sqrt{B_1^2+5{\pi }^2/4}},\forall \left(\mathbf{u},a\right)\in \overline{G}\times \left[0,\pi \right],$$

$$t_z^2\left(\mathbf{u},a\right)\ge \frac{B_1^2}{B_1^2+5{\pi }^2/4},\forall \left(\mathbf{u},a\right)\in \overline{G}\times \left[0,\pi \right],$$

$$t_z^2\left(\mathbf{u},a\right)-t_z^2\left({\mathbf{v}}_a,a\right)\ge 0,\forall \left(\mathbf{u},a\right)\in \overline{G}\times \left[0,\pi \right].$$

$$v\left(\mathbf{u},a\right)=t_z^2\left(\mathbf{u},a\right).$$

$$t\left(x,y,z,a\right)=-\underset{z}{\overset{B_2}{\int }}\sqrt{v\left(x,y,s,a\right)}ds+p\left(x,y,B_2,a\right).$$

$$v\left(x,y,z,a\right)+{\left(-\underset{z}{\overset{B_2}{\int }}\frac{v_x}{2\sqrt{v}}\left(x,y,s,a\right)ds+p_x\left(x,y,B_2,a\right)\right)}^2$$

Full text

Key Words. global convergence, semi-finite dimensional mathematical model, Carleman Weight Function, weighted Tikhonov-like functional

AMS subject classification. 35R25, 35R30

1 Introduction

We call a numerical method for a nonlinear inverse problem globally convergent if there exists a theorem claiming that this method delivers points in a sufficiently small neighborhood of the correct solution without any advanced knowledge of this neighborhood. In other words, a good first guess is not needed. We construct here a globally convergent numerical method for the Travel Time Tomography Problem (TTTP) in the 3D isotropic case with formally determined incomplete data. The TTTP is also called sometimes Inverse Kinematic Problem of Seismic, see chapter 3 in [29]. Previous publications about the 3D TTTP work only with the over determined data, in which case the number of free variables in the data $m=4$ exceeds the number $n=3$ of free variables in the unknown coefficient, see, e.g. [27, 29, 33] for the theoretical results and [38] for a numerical method. On the other hand, our data are formally determined with $m=n=3$. In the 2D case, the TTTP is always formally determined with $m=n=2,$ see [5, 26, 28] for the theory and [32] for an algorithm.

The TTTP was first considered by Herglotz [8] and Wiechert and Zoeppritz [37] in 1905 and 1907 respectively in the 1D case, due to important applications in Geophysics, see chapter 3 of [29] for some details. Recently it was discovered that, in addition to Geophysics, the TTTP also arises in the inverse problem of the recovery of the spatially distributed dielectric constant coefficient from the scattering data without the phase information in both Helmholtz equation [17] and Maxwell Equations [30].

The first globally convergent numerical method for the 3D TTTP with formally determined data was constructed in the recent work of the author [25]. In addition, the data in [25] are incomplete, so as in the current paper, see Figure 1 in [25]. The globally convergent numerical method of this paper differs from the one of [25] in one important aspect. In [25] finite differences with respect to two $x,y$ out of three $x,y,z$ spatial variables were considered. On the other hand, we consider here truncated trigonometric Fourier series with respect to $x,y$. The replacement of those finite differences with the truncated trigonometric Fourier series significantly simplifies and shortens the presentation.

In numerical methods for inverse problems, it is quite acceptable to use approximate mathematical models, see, e.g. [3, 7, 9, 10, 11]. We work below only within the framework of such a model. Our approximate mathematical model consists of two assumptions. First, we assume that a certain function associated with the solution of the governing eikonal equation can be represented via a truncated Fourier series with respect to a special orthonormal basis in $L_{2}\left(0,\pi\right)$. This basis was recently proposed by the author in [18]. Functions of that basis depend only on the position of the point source. Second, we assume, as stated above, that each component of that Fourier series can be represented via a truncated trigonometric Fourier series with respect to $x,y.$ A different approximate mathematical model for the TTTP is in [25], where finite differences instead of these truncated trigonometric Fourier series are used.

Thus, we work here with a semi-finite dimensional approximation of the original inverse problem. Assumptions of this sort are often used in numerical methods for inverse problems. Furthermore, convergence when the number of terms of truncated series tends to infinity is usually not proven, see, e.g. [7, 9, 10, 11]. Indeed, it is well known that such convergence results are very challenging to prove due to the ill-posed nature of inverse problems.

A conventional Tikhonov least squares cost functional for a nonlinear inverse problem is non convex. Hence, there is no guarantee that this functional does not have local minima and ravines, see, e.g. [31] for a good numerical example of local minima even for a rather simple inverse problem. Local minima might lead to incorrect solutions since any gradient-like method of the optimization of that functional can stop at any local minimum. To avoid the local minima, it was proposed in [4, 12, 13, 16] to construct globally strictly convex cost functionals for Coefficient Inverse Problems (CIPs) and in [15] for ill-posed problems for quasilinear PDEs. This procedure is called convexification. Next, it was established in [1] that the minimizer of such a functional exists, is unique and minimizers converge to the exact solution of the original problem as long as the noise in the data tends to zero. Finally, it was also established in [1] that the gradient projection method converges to the exact solution if starting from an arbitrary point of a certain ball in the Hilbert space of an arbitrary radius $R>0.$ Since smallness conditions are not imposed on $R$, then this is the global convergence by our above definition. After the theory was cleared in [1], a number of works were published in which numerical studies of the convexification were conducted [19, 20, 22, 23, 21]. We also refer here to the publication [2] for a different version of the convexification for a hyperbolic coefficient inverse problem with a non vanishing initial condition.

The main ingredient of the convexification is the presence of the Carleman Weight Function(CWF) in the weighted Tikhonov-like functional. CWF is the function which is involved as the weight function in the Carleman estimate for a corresponding PDE operator. Unlike this, in the current paper, so as in [25], we apply the CWF for a Volterra linear integral operator. As a result, we construct here a globally strictly convex weighted Tikhonov-like functional for the TTTP. Note that prior [25] CWFs were used only for PDE operators. Also, [25] is the first work in which a Carleman estimate is applied to the TTTP.

In section 2 we pose the inverse problem. In section 3 we present some preliminary considerations. In particular, we provide more insight than before in our orthonormal basis of [18]. In section 4 we formulate our semi-finite dimensional approximation for a certain function we work with and also describe our approximate mathematical model. In section 5 we prove Lipschitz stability result for our approximate mathematical model. In section 6 we describe our numerical method and formulate corresponding theorems. We prove these theorems in section 7.

2 Statement of the Problem

Below $\mathbf{u}=\left(x,y,z\right)$ denotes points in $\mathbb{R}^{3}.$ Let $B_{1},B_{2}>0$ be two numbers, $B_{2}>B_{1}$. Define the domain $G\subset\mathbb{R}^{3}$ as

[TABLE]

The lower and upper boundaries of $G$ are denoted as $\Gamma_{l}$ and $\Gamma_{up}$ respectively, also, $\partial_{1}G$ is the vertical boundary of $G,$

[TABLE]

[TABLE]

Let $b\left(\mathbf{u}\right)$ be the speed of sound in $\mathbb{R}^{3}$ and let $c\left(\mathbf{u}\right)=1/b^{2}\left(\mathbf{u}\right).$ Then $n\left(\mathbf{u}\right)=\sqrt{c\left(\mathbf{u}\right)}$ is the refractive index. Let $c_{0}>0$ be a number. We assume that the following conditions hold true:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We remark that the monotonicity condition (2.6) goes along well with Geophysics, see chapter 3 in [29] and [36]. Consider the Riemannian metric

[TABLE]

For two points $\mathbf{u}$,$\mathbf{u}_{0}\in\mathbb{R}^{3}$ consider the geodesic line $\Phi\left(\mathbf{u},\mathbf{u}_{0}\right)$ connecting them. Then the travel time of the sound from $\mathbf{u}_{0}$ to $\mathbf{u}$ is

[TABLE]

For each $\mathbf{u}_{0}$ the function $t\left(\mathbf{u},\mathbf{u}_{0}\right)$ satisfies the eikonal equation,

[TABLE]

[TABLE]

We assume below that the source $\mathbf{u}_{0}$ runs along an interval $I$ of a straight line,

[TABLE]

Therefore, we use a new notation below for $t\left(\mathbf{u},\mathbf{u}_{0}\right)$ as $t\left(\mathbf{u},a\right),$ where $a\in\left[0,\pi\right]$ is the parameter in (2.10).

We assume everywhere below the validity of the condition of the regularity of geodesic lines, which is used in many previous publications about TTTP.

Regularity Condition. For every pair of points $\left(\mathbf{u,u}_{a}\right)\in\overline{G}\times I$ there exists a single geodesic line $\Phi\left(\mathbf{u},\mathbf{u}_{a}\right)$ *connecting them. This line intersects the boundary $\partial G$ exactly twice: at a point * $s_{l}\left(\mathbf{u},\mathbf{u}_{a}\right)\in\Gamma_{l}$ and at another point $\widehat{s}\left(\mathbf{u},\mathbf{u}_{a}\right)\in\partial G\diagdown\overline{\Gamma}_{l}$. After reaching the point $\widehat{s}\left(\mathbf{u},\mathbf{u}_{a}\right),$ this line leaves the domain $G$ and never comes back to this domain. Also, the function $t\left(\mathbf{u},a\right)\in C^{2}\left(\overline{G}\times\left[0,\pi\right]\right)\cap C^{1}\left(\mathbb{R}^{3}\times\left[0,\pi\right]\right).$

Travel Time Tomography Problem (TTTP). Assume that conditions (2.4)-(2.7) hold. Determine the function $c\left(\mathbf{u}\right)$ for $\mathbf{u}\in G$ assuming that the following function $p\left(\mathbf{u},a\right)$ is known:

[TABLE]

Remark 2.1. *It follows from (2.8) that any geodesic line $\Phi\left(\mathbf{u},\mathbf{u}_{a}\right)$, which is originated at a point $\mathbf{u}_{a}\in I,$ is a straight line, as long as $z\in\left(0,B_{1}\right).$ Hence, the function *$t\left(\mathbf{u},a\right)$ is known for $\mathbf{u}\in\Gamma_{l},$

[TABLE]

3 Preliminaries

In this section we formulate/reformulate some results of [25] which we need in the follow up sections.

3.1 A special orthonormal basis

This basis was first introduced in [18]. Our numerical experience shows that this basis works well numerically [20, 23, 24]. Let $\beta>0$ be a number. Consider the set of functions $\left\{\psi_{k}\left(a\right)\right\}_{n=0}^{\infty}=\left\{\left(a+\beta\right)^{n}e^{a}\right\}_{n=0}^{\infty}.$ Orthonormalize these functions in the space $L_{2}\left(0,\pi\right)$ using the Gram-Schmidt orthonormalization procedure. We obtain the orthonormal basis

[TABLE]

in $L_{2}\left(0,\pi\right).$ Here $P_{n}\left(a\right)$ is a polynomial of the degree $n$. Note that we have done this orthonormalization numerically in [20, 23, 24] for three different inverse problems. It works fine up to first 15 functions $\left\{\varphi_{n}\left(a\right)\right\}_{n=0}^{14}$ in [24]. Let $\left(\cdot,\cdot\right)$ be the scalar product in $L_{2}\left(0,\pi\right).$ Consider the numbers $d_{m,n}=\left(\varphi_{m}^{\prime},\varphi_{n}\right).$ For an integer $N\geq 1$ consider the $N\times N$ matrix $A_{N}=\left(d_{m,n}\right)_{m,n=0}^{N-1}.$ Then this matrix is invertible [18].

We now present some details which were not discussed previously. Consider $Span\left(\varphi_{0}\left(a\right),...,\varphi_{N-1}\left(a\right)\right).$ This is an $N-$dimensional subspace $L_{2,N}\left(0,\pi\right)$ of the space $L_{2}\left(0,\pi\right).$ Let $D_{N}:L_{2}\left(0,\pi\right)\rightarrow L_{2,N}\left(0,\pi\right)$ be the orthogonal projection operator. Then

[TABLE]

Hence,

[TABLE]

It follows from (3.1)-(3.3) that

[TABLE]

Multiply both sides of (3.3) sequentially by the functions $\varphi_{n}\left(a\right)$ and integrate with respect to $a\in\left(0,\pi\right).$ Denote

[TABLE]

We obtain $A_{N}\left(z_{N}\right)=y^{N}.$ Hence,

[TABLE]

It follows from (3.2)-(3.5) that every function $\widetilde{f}\in L_{2,N}\left(0,\pi\right)$ can be uniquely determined from its first derivative without a knowledge of any initial condition $\widetilde{f}\left(a_{0}\right)$ for any point $a_{0}\in\left(0,\pi\right).$ In fact, this is the reason why this basis was originally introduced in [18].

However, even though the function $D_{N}f$ is sufficiently close to the function $f$ in the $L_{2}\left(0,\pi\right)-$norm for sufficiently large values of $N$, this does not imply that functions $d\left(\left(D_{N}f\right)\right)/da$ and $D_{N}\left(f^{\prime}\right)$ are close to each other in the $L_{2}\left(0,\pi\right)-$norm for $f\in C^{1}\left[0,\pi\right]$. Furthermore,

[TABLE]

Therefore, both here and in [18, 20, 23, 24, KTTP], when we recover functions of $L_{2,N}\left(0,\pi\right)$ from their derivatives without knowledge of initial conditions, we work only in the approximate sense of the orthogonal projection operator $D_{N}:L_{2}\left(0,\pi\right)\rightarrow L_{2,N}\left(0,\pi\right)$.

3.2 Estimates of functions $t_{z},t_{z}^{2}$ from the below

Let $\left(\mathbf{u},a\right)\in G\times\left(0,\pi\right)$ be an arbitrary pair of points and let $\Phi\left(\mathbf{u},a\right)$ be the geodesic line connecting points $\mathbf{u}_{a}$ and $\mathbf{u}$. Let $\Phi_{0}\left(\mathbf{u},a\right)$ be the part of this line located between the planes $\left\{z=0\right\}$ and $\left\{z=B_{1}\right\}$ and let $\mathbf{v}_{a}=\left(x_{a},y_{a},B_{1}\right)$ be the point of the intersection of $\Phi_{0}\left(\mathbf{u},a\right)$ with $\left\{z=B_{1}\right\}.$ Then $0<x_{a},y_{a}<\pi$ and $\Phi_{0}\left(\mathbf{u},a\right)$ is the segment of the straight line connecting points $\left(a,\pi/2,0\right)$ and $\left(x_{a},y_{a},B_{1}\right).$ It was shown in the proof of Lemma 4.1 of [25] that such lines can be parametrized via the variable $z$. Hence,

[TABLE]

where the number $Z=Z\left(\mathbf{u},a\right)>0$ depends on $\mathbf{u}$ and $a$ and the point $\mathbf{u}=\left(x\left(Z,a\right),y\left(Z,a\right),Z\right)\in\partial G\diagdown\Gamma_{l}.$ It was also shown in that proof that

[TABLE]

By (2.12)

[TABLE]

Lemma 3.1. Let conditions (2.4)-(2.7) be in place. *Then *

[TABLE]

* *

[TABLE]

[TABLE]

Here, (3.8) is proven in Lemma 4.1 of [25] (with slightly different notations), (3.9) follows from (3.8). In addition, (3.9) follows from (2.6), (3.6) and (3.7). Finally, (3.10) follows from (2.6) and (3.6).

3.3 Eikonal equation in an integro differential form

Estimates (3.8) and (3.9) enable us to introduce the key change of variables,

[TABLE]

By (3.8) and (3.11) $t_{z}\left(\mathbf{u},a\right)=\sqrt{v\left(\mathbf{u},a\right)}$ . Hence, using the data (2.11), we obtain

[TABLE]

Differentiating (3.12), we find formulas for $t_{x}\left(x,y,z,a\right)$ and $t_{y}\left(x,y,z,a\right).$ Substituting these as well as (3.11) in (2.9), we obtain

[TABLE]

[TABLE]

where $\left(x,y,z\right)\in G,a\in\left(0,\pi\right).$ Just as in the first step of the Bukhgeim-Klibanov method [6, 14], differentiate both sides of equation (3.13) with respect to the parameter $a$ from which the right hand side of (3.13) does not depend. Since $\partial_{a}c\left(\mathbf{u}\right)\equiv 0,$ then we obtain an integro differential equation with which we work below,

[TABLE]

[TABLE]

4 Semi-finite dimensional approximation for equation (3.14)

4.1 Convenient form of the function $v\left(\mathbf{u},a\right)$

By (3.6), (3.7), (3.10) and (3.11), we seek the solution of equation (3.14) in the following form:

[TABLE]

[TABLE]

[TABLE]

We need to obtain zero Dirichlet boundary condition at vertical sides of $G$, i.e. at $\partial_{1}G,$ see (2.3). Hence, consider the data $p\left(\mathbf{u},a\right)$ in (2.11). First, we apply (4.1)-(4.3) to get:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Based on (4.4)-(4.7), we assume that there exists a function $q\left(\mathbf{u},a\right)$ such that

[TABLE]

[TABLE]

[TABLE]

We also assume that the function $q\left(\mathbf{u},a\right)$ is known. Define the function $V\left(\mathbf{u},a\right)\in C^{1}\left(\overline{G}\times\left[0,\pi\right]\right)$ as

[TABLE]

Then (4.3), (4.5), (4.7) and (4.11) imply that the the function $V\left(\mathbf{u},a\right)$ satisfies the following boundary conditions:

[TABLE]

[TABLE]

Hence, by (4.1) and (4.11)

[TABLE]

In the right hand side of (4.14), functions $v\left(\mathbf{v}_{a},a\right)>0$ and $q\left(\mathbf{u},a\right)$ are known, the function $V\left(\mathbf{u},a\right)$ is unknown and it satisfies boundary conditions (4.12), (4.13). We focus below on the search of the function $V\left(\mathbf{u},a\right).$

4.2 Approximate mathematical model

First, we represent the function $V\left(\mathbf{u},a\right)$ via a truncated Fourier series with respect to the orthonormal basis $\left\{\varphi_{k}\left(a\right)\right\}_{k=0}^{\infty}$ of section 3.1. More precisely we assume that for $\left(\mathbf{u},a\right)\in\overline{G}\times\left[0,\pi\right]$

[TABLE]

where the coefficients $V_{n}\left(\mathbf{u}\right)$ are unknown. Next, using (4.12), we assume that for $\left(x,y,z\right)\in\left[0,\pi\right]^{2}\times\left[B_{1},B_{2}\right],n\in\left[0,N-1\right]$

[TABLE]

[TABLE]

[TABLE]

where functions $w_{n,km}\left(z\right)$ are unknown. Boundary condition (4.18) is generated by (4.13). Hence (4.14) and (4.16) imply that

[TABLE]

Approximate Mathematical Model. This model consists of the assumption that the function $V\left(\mathbf{u},a\right)$ in (4.12)-(4.14) can be represented via (4.15), (4.19) with conditions (4.17), (4.18) and the substitution of so the obtained function $v\left(\mathbf{u},a\right)$ of (4.14) in the left hand side of equation (3.14) provides zero in its right hand side. Furthermore, we assume that, in the case of noiseless data, the function $c^{\ast}\left(\mathbf{u}\right)$ obtained in an obvious manner from (4.1)-(4.19) and (3.13) is independent on the parameter $a$, satisfies the first condition (2.4) as well as conditions (2.5), (2.6) in the domain $G$ and can be extended in the entire space $\mathbb{R}^{3}$ in such a way that it will satisfy the second condition (2.4) and condition (2.7), and, also, regularity condition holds for $c^{\ast}\left(\mathbf{u}\right)$.

Remark 4.1. The suitable values of numbers $N$ and $K$ should be chosen numerically, see, e.g. [7, 10, 11, 20, 22, 23, 21, 24] for such choices for a variety of inverse problems.

Substitute (4.19) in (3.14). Then we obtain the following equation for $\left(\mathbf{u},a\right)\in G\times\left(0,\pi\right)$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Multiply both sides of equation (4.20) by $\sin\left(kx\right)\sin\left(my\right),k,m=1,...,K$ and integrate with respect to $x,y\in\left(0,\pi\right)^{2}.$ We obtain

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Denote

[TABLE]

Multiply both sides of (4.21) sequentially by functions $\varphi_{0}\left(a\right),...,\varphi_{N-1}\left(a\right)$ and integrate then with respect to $a\in\left(0,\pi\right).$ Then we obtain a system of $N$ coupled nonlinear integro differential equations with respect to the matrix $Q\left(z\right).$ The left hand side of this system is $A_{N}\left(Q_{km}\left(z\right)\right),$ where the invertible matrix $A_{N}$ was introduced in section 3.1. Multiplying by $A_{N}^{-1}$ from the left and varying $k,m=1,...,K,$ we obtain another system of coupled nonlinear integro differential equations in the form:

[TABLE]

[TABLE]

Boundary condition (4.24) follows from (4.18). Here the vector functions $F,g$ depend on $N,K$ and are twice continuously differentiable with respect to their variables as long as

[TABLE]

[TABLE]

The vector $\widehat{v}_{0}$ and the vector function $\widehat{q}\left(z\right)\in C\left[B_{1},B_{2}\right]$ linearly depend on functions $v\left(\mathbf{v}_{a},a\right)$ in (3.7) and $q\left(\mathbf{u},a\right)$ in (4.8) respectively. The vector $\widehat{p}$ depends on functions $p_{x}\left(x,y,B_{2},a\right),p_{y}\left(x,y,B_{2},a\right).$

Remarks 4.2:

1*.* With respect to (4.25) recall that by (4.9) $q\left(\mathbf{u},a\right)\geq 0$ in $\overline{G}\times\left[0,\pi\right].$ The number $d$ in (4.26) is the number of ones choice. In principle, (3.10), (4.1), (4.14) and (4.16) seem to imply that one should choose $d=0$ and then replace “$>$” in (4.25) with “$\geq$*”. However, since the function *$q\left(\mathbf{u},a\right)$ is generated by the data $p\left(\mathbf{u},a\right)$, see *(2.11), (4.4) and (4.8) and since these data might contain noise, then the noisy *$q\left(\mathbf{u},a\right)$ might be non-positive at some points. Hence, the choice (4.26) provides more flexibility in terms of noise. It follows from (3.7), (3.10) and (4.1) that (4.25) and (4.26) guarantee that $v\left(\mathbf{u},a\right)>B_{1}^{2}/\left(B_{1}^{2}+5\pi^{2}/4\right)-d=const.>0$ for $\left(\mathbf{u},a\right)\in\overline{G}\times\left[0,\pi\right].$ And we obviously need this inequality in (4.21).

2. The requirement $v\left(\mathbf{u},a\right)\geq const.>0$ is a technical condition. By the numerical experience of the author in other inverse problems, most likely, this inequality will always be satisfied in computations if using the numerical method of this paper.

3. Due to the boundary condition (4.24), the problem (4.23), (4.24) cannot be solved just as a system of nonlinear coupled Volterra integral equations.

Let $q_{1}\left(\mathbf{u},a\right)$ and $q_{2}\left(\mathbf{u},a\right)$ be two functions $q\left(\mathbf{u},a\right)$ which generate vector functions $\widehat{q}_{1}\left(z\right)$ and $\widehat{q}_{2}\left(z\right)$ respectively. Then the following estimate is easy to prove

[TABLE]

where the number $C_{1}=C_{1}\left(G,N,K\right)>0$ depends only on listed parameters.

Here and everywhere below the norm of a vector function in a conventional Banah space is defined as the square root of the sum of squares of norms of its components in that space. Thus, we do not introduce below special notations for norms of vector functions, for brevity.

Below the $NK^{2}-$dimensional vector functions $Q\left(z\right)$ have the form (4.22). Let $R>0$ be an arbitrary number and the number $d\in\left(0,B_{1}^{2}/\left(B_{1}^{2}+5\pi^{2}/4\right)\right)$ be the one chosen in (4.25) (the second Remark 4.2). Denote

[TABLE]

[TABLE]

Obviously, the set $B\left(R,q,d\right)$ is convex. Lemmata 4.1 and 4.2 can be proven similarly with proofs of Lemmata 5.2 and 5.3 of [25]* *respectively using the multidimensional analog of Taylor formula [35]. Hence, we omit these proofs. Recall that by the embedding theorem $H_{0}^{1}\left(B_{1},B_{2}\right)\subset C\left[B_{1},B_{2}\right]$ implying $B\left(R,q,d\right)\subset C\left[B_{1},B_{2}\right]$ and

[TABLE]

where the constant $C=C\left(N,K,B_{1},B_{2}\right)>0$ depends only on listed parameters.

Lemma 4.1. Let $p_{1}\left(x,y,B_{2},a\right),p_{2}\left(x,y,B_{2},a\right)\in C^{1}\left(\overline{\Gamma}_{u}\times\left[0,\pi\right]\right)$ be two functions of the data on $\overline{\Gamma}_{u}$ in (2.11). Let** $q_{1}\left(\mathbf{u},a\right)$,$q_{2}\left(\mathbf{u},a\right)\in C^{1}\left(\overline{G}\times\left[0,\pi\right]\right)$ be two functions $q\left(\mathbf{u},a\right).$ For $i=1,2$, let

[TABLE]

*where $A=const.>0.$ Let $\widehat{p}_{1}$ and $\widehat{p}_{2}$ be corresponding vectors $\widehat{p}$ generated by $p_{1}\left(x,y,B_{2},a\right)$ and $p_{2}\left(x,y,B_{2},a\right)$ respectively. Denote $\widetilde{p}=\widehat{p}_{1}-\widehat{p}_{2}.$ And similarly for $\widehat{q}_{1},\widehat{q}_{2}$ denote $\widetilde{q}=\widehat{q}_{1}-\widehat{q}_{2}.$ Then the following analog of Taylor formula is valid for any pair of vector functions $Q_{1}\in B\left(R,q_{1},d\right),Q_{2}\in B\left(R,q_{2},d\right)$ with *$\widetilde{Q}=Q_{2}-Q_{1}:$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where vector functions $F_{j},j=1,2,3,4$ depend on variables listed in them as on parameters, they are continuous with respect to these parameters as long as for $i=1,2$ the vector functions $Q_{i}\in B\left(R,q_{i},d\right)$. The function $F_{4}$ depends nonlinearly on $\widetilde{Q}.$ Furthermore, the following estimates are valid:

[TABLE]

[TABLE]

[TABLE]

*Here and below $C_{2}=C_{2}\left(G,N,K,A,R,d\right)>0$ denotes different numbers depending only on listed parameters. *

Lemma 4.2. Let conditions of Lemma 4.1 hold. Then

[TABLE]

[TABLE]

[TABLE]

where vector functions $S_{k},k=1,2,3$ have the same properties as those of vector functions $F_{k},k=1,2,3$ of Lemma 4.1.

4.3 A sufficient condition ensuring (4.25)

Even though (4.25) is a technical requirement (second Remark 4.2), we provide in this section a sufficient condition which ensures (4.25) for $d=0$. Since our numerical method is seeking the vector function $Q\left(z\right),$ then it is convenient to formulate this condition in terms of components $w_{n,km}\left(z\right)$ of this vector. The proof of Lemma 4.3 is similar with the proof of Lemma 3.1 of [25]. There is an important difference, however: functions $\sin\left(kx\right)\sin\left(my\right)$ where not used in [25].

Lemma 4.3. *Assume that for every point *$\mathbf{u}\in\overline{G}$ the function $q\left(\mathbf{u},a\right)$ *in (4.8) belongs to the subspace $L_{2,N}\left(0,\pi\right)$ (section 3.1). Let *$V\left(\mathbf{u},a\right)$ be the function defined in (4.11). Denote $q_{n}\left(\mathbf{u}\right)=\left(q\left(\mathbf{u},a\right),\varphi_{n}\left(a\right)\right)$ *and *$V_{n}\left(\mathbf{u}\right)=\left(V\left(\mathbf{u},a\right),\varphi_{n}\left(a\right)\right)$ where $n=0,...,N-1.$ Assume that condition (4.16) holds. Define vector functions $\psi^{N}\left(a\right)$ and $\varphi^{N}\left(a\right)$ as $\psi^{N}\left(a\right)=\left(\psi_{0},...,\psi_{N-1}\right)^{T}\left(a\right)$ and $\varphi^{N}\left(a\right)=\left(\varphi_{0},...,\varphi_{N-1}\right)^{T}\left(a\right),$ where functions $\psi_{n}\left(a\right)=\left(a+\beta\right)^{n}e^{a}$ were defined in section 3.1 and functions $\varphi_{n}\left(a\right)$ are obtained from $\psi^{N}\left(a\right)$ via the Gram-Schmidt orthonormalization procedure, see (3.1). Let $Y_{N}$ be the matrix of this procedure, $Y_{N}\left(\psi^{N}\right)=\varphi^{N}.$ Consider the vector functions $s\left(\mathbf{u}\right)=\left(q_{0}+V_{0},q_{1}+V_{1},...,q_{N-1}+V_{N-1}\right)^{T}\left(\mathbf{u}\right)$ and $r\left(\mathbf{u}\right)=Y_{N}^{T}\left(s\left(\mathbf{u}\right)\right),$ where

$r\left(\mathbf{u}\right)=\left(r_{0},...,r_{N-1}\right)^{T}\left(\mathbf{u}\right).$ Assume that $r_{j}\left(\mathbf{u}\right)\geq 0$ for all $\mathbf{u}\in\overline{G}$ and all $j=0,...,N-1.$ Then

[TABLE]

Proof. Let the raw number $n$ of the matrix $Y_{N}$ be $\left(k_{n,0},k_{n,1},...,k_{n,N-1}\right).$ Then

[TABLE]

Hence,

[TABLE]

[TABLE]

[TABLE]

Since $\psi_{j}\left(a\right)>0$ for all $a\in\left[0,\pi\right],$ then (4.32) follows. $\square$

5 Lipschitz Stability Estimate

Theorem 5.1. Let $R>0$ be an arbitrary number. For $i=1,2,$ *let $p^{i}\left(x,y,a\right)=p_{i}\left(x,y,B_{2},a\right)\in C^{1}\left(\overline{\Gamma}_{u}\right)\times C\left[0,\pi\right]$ and *$q_{i}\left(\mathbf{u},a\right)\in C^{1}\left(\overline{G}\right)\times C\left[0,\pi\right]$ *be two pairs of functions generated by the data $p\left(\mathbf{u},a\right)$ in (2.11) and let estimates (4.31) hold. Assume that for each of these pairs there exists a solution *$Q_{i}\left(z\right)\in B\left(R,q_{i},d\right)$ *of the problem (4.23), (4.24). Let the function *$c_{i}\left(\mathbf{u}\right)$ be the corresponding right hand side of (3.13). Then the following Lipschitz stability estimates are valid:

[TABLE]

[TABLE]

where the constant $C_{2}=C_{2}\left(G,N,K,A,R,d\right)>0$ was introduced in Lemma 4.1 and it depends only on listed parameters.

* ***Proof. **Using (4.23) and Lemma 4.2, we obtain for $z\in\left[B_{1},B_{2}\right]$ a system of linear Volterra integral equations,

[TABLE]

[TABLE]

Hence,

[TABLE]

Hence, (5.1) follows from Gronwall’s inequality. Finally, (5.2) follows from (3.13) and (5.1). $\square$

Remark 5.1. Uniqueness of the reconstruction of the function $c\left(\mathbf{u}\right)$ follows immediately from Lemma 5.1. In the proof of Theorem 5.1 we have not used the boundary condition (4.24) $Q\left(B_{1}\right)=0$ and thus came up with the system (5.3) of Volterra linear integral equations. This is because Theorem 5.1 is about the Lipschitz stability estimate rather than about a numerical method. However, when constructing solution of equation (4.23), we need to use condition (4.24). And this is done in the next section.

6 Numerical Method

In this section we construct a globally convergent numerical method which solved problem (4.23), (4.24). Since the function $c\left(\mathbf{u}\right)$ can be straightforwardly computed from the function $Q\left(z\right)$ if using (3.13), (4.1)-(4.19), then this method also reconstructs the target function $c\left(\mathbf{u}\right).$

Lemma 6.1 (Carleman estimate for the Volterra linear integral operator) [25]. *Let $\lambda>0$ be a parameter. Then the following Carleman estimate is valid with the Carleman Weight  Function *$e^{2\lambda z}:$

[TABLE]

Let the pair of functions $\left(p\left(x,y,B_{2},a\right),q\left(\mathbf{u},a\right)\right)$ satisfies conditions (4.4), (4.8), (4.10) as well as condition (4.31) in which $p_{1}=p_{2}=p$ and $q_{1}=q_{2}=q$. Let $d$ be an arbitrary number satisfying (4.26) and $R>0$ be an arbitrary number. Let $B\left(R,q,d\right)\subset H_{0}^{1}\left(B_{1},B_{2}\right)$ be the convex set defined in (4.29), also see (4.28) and (4.30). Thus, in particular, $B\left(R,q,d\right)$ consists of $NK^{2}-$dimensional functions $Q\left(z\right)\in H_{0}^{1}\left(B_{1},B_{2}\right).$

We numerically solve the problem (4.23), (4.24) via the minimization on the set $B\left(R,q,d\right)$ of the following Tikhonov-like functional with the CWF $e^{2\lambda z}$ in it$:$

[TABLE]

[TABLE]

Here $\alpha\in\left(0,1\right)$ is the regularization parameter. We use the multiplier $e^{-2\lambda B_{1}}$ to balance two terms in the right hand side of (6.1). Theorem 6.1 is the main result of this paper.

Theorem 6.1. The functional (6.1) has the Frechét derivative $J_{\lambda,\alpha}^{\prime}\left(Q\right)$ at any point of the set $B\left(2R,q,d\right).$ Furthermore, this derivative satisfies the Lipschitz continuity condition

[TABLE]

*for all *$Q_{1},Q_{2}\in B\left(2R,q,d\right),$ *where the constant $M>0$ is independent on *$Q_{1},Q_{2}.$ Most importantly, there exists a sufficiently large number

[TABLE]

* depending only on listed parameters such that the functional $J_{\lambda,\alpha}\left(Q\right)$ is strictly convex on the set $\overline{B\left(R,q,d\right)}$ for any value of $\lambda\geq\overline{\lambda}.$ In other words, the following inequality is valid*

[TABLE]

[TABLE]

[TABLE]

Remark 6.1. Although Theorem 6.1, so as other theorems of this section, is valid only for sufficiently large values of the parameter $\lambda,$ our past computational experience with the convexification demonstrates that usually once can select quite reasonable values of $\lambda\in\left[1,3\right]$ in the numerical practice [1, 19, 20, 22, 23, 21].

Our next theorem is about the existence and uniqueness of the minimizer of the functional $J_{\lambda,\alpha}\left(Q\right)$ on the set $\overline{B\left(R,q,d\right)}.$

Theorem 6.2. Let $\overline{\lambda}>1$ be the number of Theorem 6.1. Then for any $\lambda\geq\overline{\lambda}$ there exists unique minimizer $Q_{\min}\in\overline{B\left(R,q,d\right)}$ of the functional $J_{\lambda,\alpha}\left(Q\right)$ on the set $\overline{B\left(R,q,d\right)}.$ In addition, the following inequality is valid:

[TABLE]

To find the minimizer $Q_{\min}$ of Theorem 6.2, we apply the gradient projection method. Let $T:H_{0}^{1}\left(B_{1},B_{2}\right)\rightarrow\overline{B\left(R,q,d\right)}$ be the orthogonal projection operator of the space $H_{0}^{1}\left(B_{1},B_{2}\right)$ on the convex set $\overline{B\left(R,q,d\right)}.$ Let $\kappa\in\left(0,1\right)$ be a number. We choose an arbitrary point $Q^{\left(0\right)}\in B\left(R,q,d\right)$ as the starting point for iterations. The gradient projection method works with the following sequence:

[TABLE]

Theorem 6.3. Let $\lambda\geq\overline{\lambda},$ where $\overline{\lambda}>1$ is the number introduced in Theorem 6.1. Let $Q_{\min}\in\overline{B\left(R,q,d\right)}$ be the minimizer of Theorem 6.2. Then there exists a number $\overline{\kappa}=\overline{\kappa}\left(G,N,K,A,R,d,\lambda\right)\in\left(0,1\right)$ depending only on listed parameters such that for every number $\kappa\in\left(0,\overline{\kappa}\right)$ there exists a number $\theta=\theta\left(\kappa\right)\in\left(0,1\right)$ such that

[TABLE]

To prove the convergence of the sequence (6.5) to the exact solution, we assume first that there exists exact solution $Q^{\ast}\in B\left(R,q^{\ast},d\right)$ of the problem (4.23), (4.24) with idealized noiseless data $\left(p^{\ast}\left(x,y,B_{2},a\right),q^{\ast}\left(\mathbf{u},a\right)\right).$ Such an assumption is a common place in the theory of ill-posed problems, see, e.g. [3, 34]. We also assume that functions $p^{\ast},q^{\ast}$ satisfy conditions (4.4)-(4.10) as well as condition (4.31) in which $p_{1}=p_{2}=p^{\ast}$ and $q_{1}=q_{2}=q^{\ast}.$ Next, let $\delta\in\left(0,d/2\right)$ be the level of noise in the data. Hence, we assume that

[TABLE]

Let $c^{\ast}\left(\mathbf{u}\right)$ be the exact coefficient $c\left(\mathbf{u}\right)$ which corresponds to the above exact data within the framework of our approximate mathematical model of section 4.2. Combining (4.1)-(4.19) with (3.13) in an obvious manner, one can construct approximations $\overline{c}^{\left(n\right)}\left(\mathbf{u},a\right)$ for the target coefficient $c^{\ast}\left(\mathbf{u}\right).$ Since functions $\overline{c}^{\left(n\right)}\left(\mathbf{u},a\right)$ might depend on the parameter $a\in\left[0,\pi\right],$ then we average these approximations with respect to $a$, thus obtaining functions $c^{\left(n\right)}\left(\mathbf{u}\right),$

[TABLE]

Theorem 6.4. *Assume that conditions of Theorem 6.3 hold. Fix an arbitrary number $\widetilde{\lambda}\geq\overline{\lambda}$ and set *$\lambda=\widetilde{\lambda}$ in the functional $J_{\lambda,\alpha}\left(Q\right).$ Also, let the regularization parameter $\alpha=\alpha\left(\delta\right)=\delta^{2}.$ Then the following convergence estimates are valid:

[TABLE]

[TABLE]

[TABLE]

Remark 6.2.* Since $Q^{\left(0\right)}$ is an arbitrary point of the set $B\left(R,q,d\right)$ and since the size of this set $R>0$ is an arbitrary fixed number, then estimates (6.7)-(6.9) mean global convergence, as defined in the first paragraph of Introduction.*

As soon as Theorem 6.1 is established, lemma 2.1 of [1] implies Theorem 6.2. Next, if Theorems 6.1 and 6.2 are true, then theorem 2.1 [1] leads to Theorem 6.3. Therefore, we are left to prove only two theorems: 6.1 and 6.4.

7 Proving Theorems 6.1 and 6.4

7.1 Theorem 6.1

The goal of this section is to prove Theorem 6.1. Let $Q_{1},Q_{2}\in B\left(R,q,d\right)$ be two arbitrary points of this set. Denote $\widetilde{Q}=Q_{2}-Q_{1}\in H_{0}^{1}\left(B_{1},B_{2}\right).$ Then $Q_{2}=Q_{1}+\widetilde{Q}$ and $\left\|\widetilde{Q}\right\|_{H^{1}\left(B_{1},B_{2}\right)}<2R.$ Hence, by (4.30)

[TABLE]

Setting in Lemma 4.1 $\widehat{p}_{1}=\widehat{p}_{2}=\widehat{p}$ and $\widehat{q}_{1}=\widehat{q}_{2}=\widehat{q}$, we obtain for $z\in\left(B_{1},B_{2}\right)$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $D_{1}$ and $D_{2}$ are continuos functions of their variables, the function $D_{2}$ depends nonlinearly on $\widetilde{Q}$ and the following two estimates are valid:

[TABLE]

[TABLE]

Denote

[TABLE]

Hence,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

First, we estimate from the below all lines of (7.5), except of the first two. Using (7.1), (7.3) and (7.4), we obtain for the third line:

[TABLE]

Next, using Cauchy-Schwarz inequality, we obtain for the fourth line:

[TABLE]

Similarly for the fifth, sixth and seventh lines,

[TABLE]

[TABLE]

[TABLE]

Denote $Z$ the sum of all lines of (7.5), except of the first two. Hence, (7.6)-(7.9) imply

[TABLE]

The second line of (7.5) is linear with respect to $\widetilde{Q}\left(z\right).$ Hence, (6.1) and (7.5) lead to

[TABLE]

where $I_{lin}\left(\widetilde{Q}\right)$ depends linearly and $I_{nonlin}\left(\widetilde{Q}\right)$ depends nonlinearly on $\widetilde{Q}.$  More precisely, the linear part is:

[TABLE]

[TABLE]

[TABLE]

where $\left[\cdot,\cdot\right]$ denotes the scalar product in the space $H^{1}\left(B_{1},B_{2}\right)$ of $NK^{2}-$dimensional functions. As to the nonlinear part, it satisfies estimates (7.5), (7.14), where (7.13) is obvious and (7.14) follows from (7.10):

[TABLE]

[TABLE]

[TABLE]

It follows from (7.2), (7.4) and (7.12) that

[TABLE]

Hence, by Riesz theorem there exists unique element $Y\in H_{0}^{1}\left(B_{1},B_{2}\right)$ such that

[TABLE]

In addition, it follows from (7.11), (7.13) and (7.15) that

[TABLE]

Hence, $I_{lin}\left(\widetilde{Q}\right)=\left[Y,\widetilde{Q}\right]$ is the Frechét derivative of the functional $J_{\lambda,\alpha}\left(Q\right)$ at the point $Q_{1},$ i.e.

[TABLE]

We omit the proof of the Lipschitz continuity property (6.2) of the Frechét derivative $J_{\lambda,\alpha}^{\prime}\left(Q\right)$ since this proof completely similar with the proof of Theorem 3.1 of [1].

It follows from (7.11), (7.14) and (7.16) that

[TABLE]

[TABLE]

[TABLE]

By Lemma 6.1

[TABLE]

Choose $\overline{\lambda}>1$ so large that $C_{2}/\left(2\lambda\right)<1/8$ for all $\lambda\geq\overline{\lambda}.$ Hence, (7.17) implies that

[TABLE]

[TABLE]

7.2 Theorem 6.4

In this section we prove Theorem 6.4. Setting in Lemma 4.2

[TABLE]

[TABLE]

we obtain

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Since

[TABLE]

then (6.6) and (7.18)-(7.20) imply that

[TABLE]

To indicate the dependence of the functional (6.1) on vector functions $p\left(x,y,B_{2},a\right)$ and $q\left(\mathbf{u},a\right),$ we temporary denote this functional as $J_{\lambda,\alpha}\left(Q,\widehat{p},\widehat{q}\right).$ Recall that $\alpha=\delta^{2}.$ Since the number $\lambda=\widetilde{\lambda}\geq\overline{\lambda}$ is not varied, then we can estimate $\exp\left[2\lambda\left(B_{2}-B_{1}\right)\right]\leq C_{2}.$ Hence, using (6.1) and (7.21), we obtain with a different constant $C_{2}>0:$

[TABLE]

Next, by (6.3)

[TABLE]

[TABLE]

Hence, using (7.22), we obtain

[TABLE]

By (6.4)

[TABLE]

Using (7.23) and (7.24) we obtain (6.7).

Next, the triangle inequality and Theorem 6.3 give:

[TABLE]

[TABLE]

[TABLE]

Thus,

[TABLE]

Combining (7.25) with (6.7), we obtain (6.8). Finally, estimate (6.9) follows (6.8) and the construction of functions $c^{\left(n\right)}\left(\mathbf{u}\right)$ described in section 6. $\square$

Acknowledgment

This work was supported by US Army Research Laboratory and US Army Research Office grant W911NF-19-1-0044.