{'fraction_non_alphanumeric': 0.031192771326799623, 'fraction_numerical': 0.011470890100823087, 'mean_word_length': 4.631459655485041, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Data-driven solutions of ill-posed inverse problems arising from doping reconstruction in semiconductors

The non-destructive estimation of doping concentrations in semiconductor devices is of paramount importance for many applications ranging from crystal growth, the recent redefinition of the 1kg to defect and inhomogeneity detection. A number of technologies (such as LBIC, EBIC and LPS) have been developed which allow the detection of doping variations via photovoltaic effects. The idea is to illuminate the sample at several positions, and detect the resulting voltage drop or current at the contacts. We model a general class of such photovoltaic technologies by ill-posed global and local inverse problems based on a drift-diffusion system which describes charge transport in a self-consistent electrical field. The doping profile is included as a parametric field. To numerically solve a physically relevant local inverse problem, we present three different data-driven approaches, based on least squares, multilayer perceptrons, and residual neural networks. Our data-driven methods reconstruct the doping profile for a given spatially varying voltage signal induced by a laser scan along the sample's surface. The methods are trained on synthetic data sets (pairs of discrete doping profiles and corresponding photovoltage signals at different illumination positions) which are generated by efficient physics-preserving finite volume solutions of the forward problem. While the linear least square method yields an average absolute error around 10%, the nonlinear networks roughly halve this error to 5%, respectively. Finally, we optimize the relevant hyperparameters and test the robustness of our approach with respect to noise.

The LPS method, proposed in [24,8] and schematically shown in Figure 1 right panel, is especially useful in the context of crystal growth, where it is virtually impossible to predict the quality (i.e., the symmetry) of a semiconductor crystal during its growth in a furnace. During the growth process, thermal fluctuations near the the solid-liquid interface introduce local fluctuations (or striations) in the doping profile. LPS detects such doping inhomogeneities non-invasively at wafer-scale and room temperature. It is especially suitable for low doping concentrations (10 12 cm −3 to 10 16 cm −3 ).
Mathematically speaking, all of the discussed technologies result into inverse problems. The forward problem assumes that we know the doping profile and a set of laser spot positions. We then want to know the corresponding (laser spot dependent) photovoltage signals at the contacts. The inverse problem, on the other hand, assumes we have measured photovoltage signals at the contacts for a laser scan across the sample. Then we want to reconstruct the doping profile at least in the probing region but ideally in the whole domain. We model a general class of photovoltaic technologies introduced above by ill-posed global and local inverse problems based on a drift-diffusion system which describes charge transport in a self-consistent electrical field, the so-called van Roosbroeck system. The doping profile is included as a parametric field. We point out that a key difficulty in the mathematical modeling and numerical simulation is that the probing area where the laser or electron beam operates is significantly smaller (even up to one spatial dimension) than the region in which we would like to recover the doping profile. This means that the solution of such incomplete inverse problems cannot be unique and we have to look for appropriate minimizers or rely on some physically meaningful assumptions on the doping regarding periodicity or variation in only one spatial dimension.
Ill-posed inverse (PDE) problems have been studied for a long time from an analytical and a numerical point of view. We refer to the general overviews [40,15] and the references therein. Burger, Markowich and others analyzed inverse semiconductor problems [21,5]. Since the analysis is challenging, they often consider linearized/unipolar settings, no recombination/generation, only small external biases, linear diffusion as well as standard Dirichlet-Neumann boundary conditions. In [5], they discuss identifiability of the doping profile from capacitance, reduced current-voltage or laser beam induced current (LBIC) data. As for the doping reconstruction, previous numerical methods relied on optimization [33], the level set [21] or the Landweber-Kaczmarz methods [5]. Especially, the latter proved to be unstable and costly. For this reason we propose in this paper data-driven approaches to solve the local inverse photovoltage PDE problem.
In particular, we focus on three data-driven approaches: First, noting that under special conditions the photovoltage signal is related to the doping profile by linear operations, namely differentiation and convolution, we try a classical least squares approach. Second, allowing also a nonlinear relationship between doping profile and photovoltage signal, we train multilayer perceptrons. Finally, we adjust more advanced ResNets to our specific setup to solve the inverse photovoltage problem. Any data-driven technique heavily relies on the amount and quality of the training data. Since, for example, growing crystals is exceptionally costly, we are not able to generate large real-life datasets. For this reason we generate physics-preserving synthetic data (measured signals and corresponding doping profiles) via a fast and efficient implementation of the forward PDE model [8] which relies on the Voronoi finite volume discretization described in [9]. As discussed above, in the context of LPS it has been shown that the forward model nicely encompasses three main physically meaningful features [8]. The flux discretization is handled by ideas of Scharfetter and Gummel [37]. Compared to an implementation based on commercial software [18,16], our open-source code reduces the simulation time of the forward model by two orders of magnitude. Finally, we will study the robustness with respect to noise and carefully tune the hyperparameters. We introduce noise in the amplitude as well as in the wave lengths and phase shifts of the doping fluctuations. While the results for noisy data are worse than for clean data, our approach appears to be relatively robust with respect to noise.
The literature on how deep neural networks maybe be used to solve PDEs has been rapidly increasing in recent years. A notable numerical approach is called physics informed neural networks (PINN); see [35] and references therein for more details. The key idea is to replace classical statistical loss functions with PDE residuals. Unfortunately, in our case this is not directly feasible since solving the forward problem requires knowledge of the doping profile in the entire three-dimensional domain. For the inverse problem, however, the doping profile may only be reconstructed within a two-or even one-dimensional subset. Other approaches include supervised deep learning algorithms [36,28] to efficiently approximate quantities of interest for PDE solutions and [25] and reference therein to accelerate existing PDE numerical schemes. In order to solve high-dimensional PDEs, we refer to [12,11]. In terms of recovering parameters in governing PDEs by solving inverse problems, we refer to [7,29,26,31,39]. In particular, [7] considers an inverse scattering problem for nano-optics. How to recover generalized boundary conditions is studied in [29]. Active learning algorithms are also proposed to approximate associated PDE constrained optimization problems [26]. In [31], the authors provide model-constrained deep learning approaches for inverse problems. Such strategies are applied in [39] to reconstruct the second-order coefficients in elliptic problems. Especially in the context of semiconductors, it is often not clear which material parameters such as life times or reference densities are actually correct. For this reason, the authors in [20] used machine learning techniques to estimate material parameters in the context of organic semiconductors.
The remainder of the paper is organized as follows: In Section 2, we introduce general (forward) PDE models for general electromagnetic source terms and doping profiles based on the van Roosbroeck system which describes charge transport in a self-consistent electrical field. In Section 3, we present global and local inverse photovoltage models for generic electromagnetic source terms. In Section 4 we present the three data-driven approaches we use to solve the local inverse photovoltage problem for a specific LPS setup. We conclude with a summary and an outlook in Section 5.
2. Forward photovoltage model. In this section, we first describe the drift-diffusion charge transport model and then the circuit model which represents the volt meter, modeled by boundary condition. When combining both models, we are able to formulate the forward photovoltage model which may be used to predict a measured signal (either a current or a voltage) for given doping profile and source of electromagnetic radiation.

The van Roosbroeck model. The semiconductor crystal is modeled by a bounded
domain Ω ⊂ R 3 in which two charge carriers evolve: electrons with negative elementary charge −q, and holes with positive elementary charge q. The doping profile is given by the difference of donor and acceptor concentrations, We assume that C is a bounded function, and we call C ⊆ L ∞ (Ω) the space of all admissible doping profiles.
We describe the charge transport within the crystal in terms of the electrostatic potential ψ(x), and quasi-Fermi potentials for electron and holes, ϕ n (x) and ϕ p (x), respectively. The current densities for electrons and holes are given by J n (x), J p (x). These variables shall satisfy the so-called van Roosbroeck model in which the first equation, a nonlinear Poisson equation is self-consistently coupled with two continuity equations [27] −∇ · (ε∇ψ) = q(p − n + C(x)), In the van Roosbroeck model (2.1), the permittivity of the medium is denoted by ε and the mobilities of electrons and holes are respectively indicated by µ n and µ p . Assuming so-called Boltzmann statistics, the relations between the quasi-Fermi potentials and the densities of electrons and holes, n and p respectively, are given by Here, we have denoted the conduction and valence band densities of states with N c and N v , the Boltzmann constant with k B and the temperature with T . Furthermore, E c and E v refer to the constant conduction and valence band-edge energies, respectively. The semiconductor is considered in equilibrium if the quasi-Fermi potentials vanish, ϕ n = ϕ p = 0. In this case only the nonlinear Poisson equation in (2.1) is solved for the equilibrium electrostatic potential ψ eq . The corresponding equilibrium charge carrier densities n eq and p eq satisfy n eq p eq = n 2 i , where n i is the so-called intrinsic carrier density, defined via the relationship The recombination term H is the sum of the direct recombination, the Auger recombination and the Shockley-Read-Hall recombination, that is respectively: Different types of crystal samples (such as silicon, germanium, gallium arsenide) have different life times τ n , τ p , and reference densities n T , p T . The electromagnetic source (a laser or an electron beam) is modeled by the generation term G(x; x 0 ). When the laser hits the crystal at the point x 0 := (x 0 , y 0 , z 0 ) T , some photons are impinged and create electron-hole pairs, resulting in a generation rate defined as follows where S(x) is the shape function of the laser (normalized by R 3 S(x)dx = 1), while κ is a constant given by κ := P λ L h (1 − r). Here, P is the laser power, λ L is the wave length of the laser, h is the Planck constant, and r is the reflectivity rate of the crystal.
We assume that area of influence of the electromagnetic source decays exponentially fast from the incident point x 0 . In particular, we take a laser profile function S defined as Here σ L is the laser spot radius, while d A is the penetration depth (or the reciprocal of the absorption coefficient), which heavily depends on the laser wave length. Figure 2 shows a typical configuration, where the laser beam hits the crystal on the top surface, and shows the exponential decay of the laser shape function S. A prototypical setting for photovoltage measurements is given by a cuboid sample attached to a voltmeter with resistance R, and the electromagnetic source aiming at its top surface.

2.2.
Boundary conditions. The PDE system (2.1) is supplemented with Dirichlet and Neumann boundary conditions. The boundary ∂Ω is the union of two disjoint parts Γ N and Γ D . On Γ N , we assign Neumann boundary conditions where ∂/∂ν = ν · ∇ is the normal derivative along the external unit normal ν. On Γ D , we assign Dirichlet-type boundary conditions. This type of boundary condition models so-called ohmic contacts [38]. We suppose that there are two ohmic contacts, i. e. Γ D = Γ D1 ∪ Γ D2 . The ohmic boundary conditions can be summarized by where ψ 0 is the local electroneutral potential which one obtains by solving the Poisson equation for ψ for a vanishing left-hand side. The terms u Di denote the contact voltages at the corresponding ohmic contacts; for the forward model their difference is a-priori unknown. We define a reference value u ref of the potential and set it to The total electric current i D flowing through the ohmic contact Γ D2 is defined by the surface integral The current i D is a function of the contact voltage u D2 , of the laser spot position x 0 , and of the doping profile C.
In order to close the system, we model the voltage/ampere meter as a simple circuit having a resistance R. This structure is visualized in Figure 1. When the external circuit is connected to the crystal, the boundary condition u D2 can no longer be chosen arbitrarily and must balance the voltage difference at the voltmeter due to the current i D , which is equal to the product R i D . A generalized theory of this coupling can be found in [2] and in the references therein.
Let U be the function space of all possible photovoltage signals as a function of the laser spot position x 0 that can be measured at the contacts for any given admissible doping profile C ∈ C. We define the laser-voltage (L-V) map as precisely the map that associates to each laser spot location x 0 and doping profile C the corresponding photovoltage signal u P which satisfies the second ohmic boundary condition in (2.6) and balances the voltage difference across R, that is The laser-voltage map (2.8) corresponds to a nonlinear version of Ohm's law at the contact and constitutes an implicit equation for the Dirichlet boundary condition u D2 of the van Roosbroeck system (2.1). An existence result for solutions of the laser-voltage map is provided in [6], with the assumption that the generation rate G is small enough. For a spatially varying doping profile C, the laser-voltage map u P (x 0 , C) may be different for each laser spot location x 0 .

The global forward photovoltage problem.
Although the presented forward PDE model is well-posed for low laser power intensity and for all laser spot positions x 0 ∈ Ω, not all positions in Ω can be illuminated by a laser. The major limiting factor is the fact that the laser power intensity decays exponentially away from the surface of the crystal sample, leaving in fact as the only feasible positions those that are near or on the surface of the crystal sample.
We denote with Σ ⊆ Ω the set of all admissible laser spot positions x 0 . Our global forward photovoltage problem then associates to each doping profile C the function u that maps each laser spot position x 0 ∈ Σ to the corresponding photovoltage signal u P (x 0 , C), i.e.,

Inverse photovoltage problems.
From an industrial perspective, more interesting than the forward photovoltage problem is the inverse one. How do doping inhomogeneities influence the measured voltage difference u P , defined in the previous section? Suppose we have measured the photovoltage signal for several different laser spot positions x 0 , how does the doping profile look like that leads to this signal? Since the doping profile enters as a volumetric source term defined in the whole domain in the van Roosbroeck system (2.1), but we can only probe part of the domain with the laser signal, answering such a question implies solving an ill-posed inverse problem: different global doping profiles may lead to the same signal. Hence, we will first formulate an idealized global inverse photovoltage problem and then a practically relevant local inverse photovoltage problem.

Global inverse photovoltage problem.
Ideally, we would like to find the global doping reconstruction operator F := U −1 : U → P(C), that, for a given photovoltage signal measurement u, returns the preimage U −1 (u). Here we indicate with P(C) the power set of C, i.e., the set of all possible subsets of C. The global inverse photovoltage problem reads that is, F is a function from U to P(C) (i.e., a multi-valued function of the signal u) such that

Local inverse photovoltage problem.
In practice, however, the construction of the operator F (even the construction of an approximation of F ) is nontrivial, and we would like to simplify our problem by building a local inverse problem that is better posed. A key point when defining a local inverse problem, is the fact that we cannot probe the entire domain Ω with the laser, but only the subset of all possible laser spot positions Σ near or on the surface of Ω. In the following subsections, we will make some assumptions on the structure of technologically feasible doping profiles. Technological feasibility depends on the growth process, on the technology used to inject doping in the crystal, and on the specific photovoltage technology. Here, we assume variation only in the x-direction. That is, C(x) = C(x). This class of doping profiles is relevant, for example, for crystal growth, where striations in the doping profile indicate fluctuations of the temperature field during the growth process, which are dominant in the growth direction.
We define the set of technologically feasible (TF) doping profiles and the set of corresponding signals as With these assumptions on the doping profile, it is reasonable to presume that the photovoltage signal will be influenced more by variations of the doping profile in the vicinity of the laser spot position, which is where most of the charge carriers are generated, and therefore we define the local inverse photovoltage problem by restricting the global inverse problem (3.1) both in terms of technologically feasible doping profiles, as well as in terms of probing domain: where the restriction operator | is applied to each element in the set. The goal of the local inverse problem is to reconstruct doping profiles only in the probing area Σ and not on the whole domain Ω. However, the underlying assumption is that the local doping reconstruction will give us information in a neighborhood of Σ, or even on all of Ω, when one makes additional assumptions (for example on doping periodicity). In general, the set C TF,u,Σ = F (u) is much smaller than the set C u := F (u), since it restricts the family of admissible doping profiles to C TF , and discards all of the information outside of Σ. Two doping profiles in C u correspond to the same element in C TF,u,Σ as soon as their restrictions on Σ coincide. In principle, it should be possible to formulate a well-posed local inverse problem by further reducing the admissible doping profiles C TF until the set C TF,u,Σ contains just one single element for any admissible input signal u ∈ U (C TF ). We do not know, however, if this procedure is feasible, and what the necessary and sufficient conditions on C TF and on Σ are to ensure that the local inverse problem F (u) is well-posed. We leave these questions for future investigations, and concentrate on numerical approximations of the local inverse problem based on well-posed finite dimensional approximations of F . In Section 4 we propose three different strategies (of increasing complexity) to build an approximate inverse operator F h starting from a collection of doping profiles restricted to a discrete set of probing points Σ h and their corresponding discrete photovoltage signals.
Remark 3.1. The numerical approximations that we construct always produce a unique answer for each finite sampling of a signal u in U TF . We do not attempt to construct the full set C TF,u,Σ of all possible doping profiles that would generate u, but only provide the sampling of a single doping profile C, which is, in some sense, a representative of the set C TF,u,Σ . The ill-posedness of the local inverse problem and the impossibility to recover the full doping profile C after having solved the approximate inverse problem is one of the reasons why we cannot use more modern techniques (such as Physics Informed Neural Networks (PINN) [35]) that would exploit the residual of the forward problem to improve the construction of an approximate inverse operator.
4. Data-driven approximation of the inverse photovoltage problem. Inverse problems for charge transport equations have often been tackled with standard techniques (see, for example, [5,21,33]). However, in other fields such as image recognition or weather prediction, data-driven approaches have gained significant momentum to solve a variety of inverse problems (see, for example, [3,22,1]).
In order to formulate a data-driven approach for the numerical approximation of the operator F , we leverage the well-posedness of the forward problem U and its discrete approximation described in [8] to formulate a discrete inverse problem F h as a well-posed minimization problem in a finite dimensional space.

Discrete local inverse problem.
Let Σ ⊆ Ω be a subdomain. We assume we sampled both the photovoltage signal u P and the doping profile C at n ∈ N discrete laser spot positions which shall be contained in the discrete mesh Σ h ⊆ Σ. With the boldface symbols u and C we indicate the discrete samplings of the signal u P (·, C) and the doping profile C evaluated on Σ h , respectively. Notice that the generation of a single pair of signal and doping samples (u, C) ∈ R n×2 from an admissible doping profile C ∈ C TF requires in fact the solution of n discrete van Roosbroeck systems (2.1), one for each laser spot position The local discrete inverse problem F h can be interpreted as a (generally nonlinear) function that maps a vector of signal measurements u ∈ R n to a vector of doping profile values C ∈ R n : Ideally, we would like to build F h in such a way that for all u in U TF , there exists a C ∈ C TF,u,Σ = F (u) such that F h (u := u| Σ h ) = C := C| Σ h . However, this procedure suffers from the non-uniqueness of F (u), which we mitigate by constructing F h through a well-posed minimization problem.
In practice, we build F h by minimizing the mean square error loss consisting of a collection of N train pairs of signal samples u j and doping samples C j . The generation process of the data is detailed in Subsection 4.3.
The quality of the approximation of F h is evaluated by measuring the prediction capabilities of F h on a test dataset that is independent of the training dataset. For the evaluation of the prediction quality, we use a perhaps more physically relevant average ∞ norm, and define the test error as We develop three approaches to build F h : • Least squares (LS): we approximate F h (u) by the matrix vector product of an n × n matrix A h with u. We find the matrix by minimizing the MSEL error defined in (4.2) over all n × n matrices; • Multilayer perceptron (MLP): we approximate F h (u) by a multilayer perceptron, and use Stochastic Gradient Descent (SGD) to minimize the MSEL; • Residual neural network (ResNet): we approximate F h (u) by a residual neural network, and use SGD to minimize the MSEL.

An industrial application: the LPS setup.
We introduce a specific LPS setup for a silicon sample, discuss how and why we generate the data from the forward model as well as solve the corresponding local inverse photovoltage problem via the three data-driven approaches introduced in Subsection 4.1.
We consider a silicon cuboid of the form Ω : , 0] with length = 3 mm, width w = 0.5 mm and height h = 5 × 10 −5 mm, ohmic contacts at x = − /2 (corresponding to Γ D1 ) and x = /2 (corresponding to as Γ D2 ), see Figure 1. We focus on a lateral photovoltage scanning (LPS) table-top setup, see Figure 1 right panel. The silicon parameters needed for the corresponding forward model from Section 2 can be found in [8]. We simulate the charge transport in the two-dimensional plane where y = 0 and reconstruct the doping profile along a line, namely the subdomain Σ = {(x, 0, 0) T ∈ Ω : x ∈ [0, k /2]} with k < , ensuring we are sufficiently far away from the boundaries x = ± /2 to reduce boundary effects. The laser scan then produces a uniform partition 0 = x 1 < . . . < x n = k /2 with mesh size ∆x. The resulting laser spot positions are given by In the following numerical simulations, we set n = 1200 and k = 0.4 mm.

Data generation.
For our data-driven approaches, we need a suitable amount of data; unfortunately, generating enough data from real life requires large budgets. Instead, we will generate synthetic data from the forward model that we have described in Section 2. This data is physically meaningful in the sense that our finite volume approximation of the LPS problem correctly incorporates physically meaningful behavior, namely i) the signal intensity depends linearly on local doping variations for low laser powers, ii) the signal saturates for higher laser intensities due to a screening effect, and iii) the signal depends logarithmically on moderate laser intensities [8].
As pointed out in Subsection 3.2, we consider a family of doping profiles that vary only along the x direction, and that are constant in both y and z directions. To generate our synthetic datasets, we need an algorithm that produces fluctuating doping profiles. Since within the LPS framework the doping profile is roughly periodic, see [8], we assume a superposition of sinusoidal functions from which we randomly sample physically reasonable amplitudes and wavelengths. Therefore, we define In our simulation, we set the number of sine terms to N b = 5, which seems to strike a good balance between complexity and real-life situations for striations in doping profiles.
To generate an element {u, C} of our dataset, we randomly sample a parameter β with which we generate a continuous function C(x, β) and we compute the discrete counterpart of u := U (C) defined in Equation (2.9) using the finite volume scheme developed in [19]. Finally, we restrict both C and u to the discrete laser spot mesh Σ h .
In realistic physical scenarios, the doping profiles contain noises, and they cannot be represented exactly by Equation (4.4). To simulate such scenario, we perturb some of the doping profiles and generate an additional noisy dataset, used to test the robustness of our networks in the presence of noise. In what follows, we say that a dataset is "noisy" or "clean" if C has or has not been perturbed. The algorithm used to generate noisy datasets is described in Appendix A.
In the following, we summarize the datasets used to construct and test the learning models introduced in Section 4.4-4.6. The datasets are available on the repository [34].
CleanDataSET. Consisting of: • CleanTrainingDataSET: a dataset with N train = 22, 500 clean samples used to train our initial model; • CleanTestDataSET: a dataset with N test = 7, 500 clean samples. With this dataset we test the performance of our models trained by CleanTrainingDataSET. Note that testing errors generated by this dataset (cf. (4.3)) implies the quality of a trained model. This will help use tune the corresponding hyperparameters; • CleanValidationDataSET: a dataset with the same size of the CleanTestDataSET that will be used to validate the performance of our model on the clean case after performing the tuning.

NoisyDataSET. Consisting of:
• NoisySmallTestDataSET: again a dataset of N test = 7, 500 samples, but this time with noise. This is used to test the robustness of the models trained by the CleanTestDataSET; • NoisyTrainingDataSET: a dataset with N train = 150, 000 noisy samples that we use to train our models to be robust to noise; • NoisyTestDataSET: a dataset with N test = 50, 000 noisy samples. This is the dataset that we use to test the performances of our models trained with the NoisyTrainingDataSET. We provide a qualitative analysis of CleanTrainingDataSET and NoisyTrainingDataSET in Figure 3, where we compute the Singular Value Decomposition (SVD) of both datasets interpreted as n × N train matrices. The number of dominant singular values is a crude indication of the actual dimension of the spaces C TF , and U TF , and shows that, roughly, their dimension is close to 50 when using the clean data generation described above, see Figure 3 left panel. The two dimensions mismatch significantly when adding noise, see Figure 3 right panel. While the dimension of the clean photovoltage signals follows roughly that of the clean doping profiles, the same cannot be said for the noisy cases. This difference is a hint that the inverse operator F is not well-posed, due to a mismatch in the dimension of the input and output spaces.

Least squares.
In [8] the authors showed that, under certain restrictive conditions, the operator F can be approximated by a linear one. Indeed, in case of an n or p doped semiconductor that varies only along the x direction, we can show that If the support of the laser source is larger than the wavelength of oscillation of the doping profile, measuring the voltage difference may actually result in a signal that does not capture a single oscillation, but a floating average. We could represent this floating average of doping fluctuations by a convolution of the doping gradient with some (unknown) profile f , i.e.
Hence, it may seem plausible to relate the photovoltage signal to the doping profile via linear operations such as convolution and differentiation. However, it is not clear whether this profile function f is actually independent from the doping fluctuations, the shape of the laser beam or the laser spot position. The least square analysis is useful to understand how the inverse problem behaves, and what type of operation the inverse operator F h (u) = A h u performs. Let us express the training data for the photovoltage signals and doping profiles generated according to Subsection 4.3 with the two dense matrices U n×Ntrain and C n×Ntrain . Then we wish to solve the least square problem (4.7) A h = arg min B∈R n×n In Figure 4 we show A h obtained with CleanDataSET (left) and with NoisyDataSET (center). We observe a highly non-local behavior of F h . In the CleanDataSET case, this behaviour is more pronounced than in the NoisyDataSET case.
If we compare the singular values of the matrix A h in Figure 4 (right) with those of the two datasets in Figure 3, we observe that the dominant singular values of the least square matrix computed with the CleanDataSET are the first thirty (similar to what happens in the singular values of the CleanDataSET itself in Figure 3), while those that are most relevant for the matrix constructed with the NoisyDataSET are the first one hundred. While this information is only qualitative, it does show that a non-negligible part of the local inverse problem can be approximated relatively well by the linear operator F h (u) = A h u, with an intrinsic dimension of around one hundred, when including noise, and much smaller when considering a clean dataset. The statistical distribution of the test error in the CleanDataSET (i.e., the error defined in (4.3) for the CleanTestDataSET) is reported in Table 2 and depicted in the left panel of Figure 5. It shows an error which is, on average (orange dashed line), around 9%. Even though such a model presents a relatively high error, we observe that it is still capable of capturing the overall qualitative behavior of the doping profiles (see Figure 6).
Next, we check how robust our networks are with respect to noise. As a first test, we use the linear model generated with the CleanTrainingDataSET to predict the doping of the samples in the NoisyTestDataSET. Since the networks have not been trained with noise, we observe a significant deterioration in the quality of the predictions, shown in the middle plot of Figure 5, where the average error is now around 20%.
We believe that the main reason why the error increases so much is the fact that the noise introduces a shift in the average of the doping function that can not be physically estimated by our setup. In other words, Equation (roughly) equal to 1 · 10 16 = C 0 on the overall domain. When we introduce noise, this assumption does not hold anymore and we have no way to predict whether the average of the doping is the one we expect in the entire crystal sample. This is also coherent with the qualitative analysis in Equation (4.5), where we show that the value of the current is mostly related to the local variation of the doping and not to its average value. Removing the average from both the output of the model, and from the reference doping during testing leads to the results in the right plot of Figure 5 (bottom), where the error drops from 20% to 15%. Indeed, we expect still a larger error w.r.t. to the histogram presented in Figure 5, since the test dataset includes noisy data which are not included in the training dataset.
A few examples of predictions of the doping profile on Σ h with noisy input data are shown in Figure 7.
An improvement of the error for the least square problem can be obtained by performing the training on the NoisyTrainingDataSET. The histogram of the error on the NoisyTest-DataSET is shown in the right plot of Figure 5, where the error is now around 10% (see Table 2 for the details). It becomes apparent that training with noisy doping profiles also significantly reduces the need to remove the average doping value.
In summary, the least square model is able to predict doping profiles with an average error of around 10% for clean test data or noisy but properly average adjusted test data. However, it is possible to further improve the results by introducing nonlinearities in our models, for example using neural networks. We focus on multilayer perceptrons and residual neural networks in the following sections.

Multilayer perceptron.
While the least square approach is a good starting point, its efficiency is only good if the inverse operator is close to a linear operator. For inverse problems associated to the van Roosbroeck system, this is not necessarily the case, and one may choose to introduce some nonlinearities in the approximation of F h . Multilayer perceptrons [10] are among the most widely used feedforward neural networks, and they are the first natural candidate for general nonlinear function approximations. We consider networks consisting of a down-sampler There are fewer than 6 4 · 9 2 configurations because some of them would lead to a negative amount of neurons. The multilayer perceptrons are implemented in PyTorch [32].
For a given configuration of the neural network, we apply the stochastic gradient descent (SGD) algorithm (without momentum and weight decay) to find the MLP that minimizes the mean square error loss defined in Equation (4.2). For each network configuration, we fix the learning rate to a value that is chosen randomly in the interval [10 −3 , 1] and the batch size to 64. The number of training epochs is set to 1, 000.
We randomly selected 10, 000 configurations where we vary both the structure of the MLP (choosing one of the 71,118 possible MLPs) as well as the learning parameters of SGD. To reduce the computational cost associated with the training of all the resulting neural networks, we apply the Asynchronous Successive Halving algorithm (ASHA) [23], a multi-armed bandit algorithm that has been optimized to run on a massive amount of parallel machines. The ASHA algorithm discards the worst 50% performers of the neural network population at each check point, and proceeds with the training only for the most promising neural networks. Our implementation is based on the Ray library [30], which also takes care of distributing the jobs and coordinating the execution among different nodes during the computation.
The multilayer perceptron with the smallest ∞ testing error defined in Equation tistical distribution of the error is reported in the middle part of Table 2, and depicted in the left plot of Figure 8. The ∞ error is, on average, around 4.67%. Even with such a simple neural network, we reduce the average error with respect to the least square model by of a factor two. Also, we obtain a much better statistical distribution, clearly shown in Figure 8. Furthermore, the resulting MLP is robust to noise, provided that we properly filter the average of the doping profile. A comparison between the statistical distribution of the errors (both with and without average) as well as the corresponding statistical distributions are shown in Table 2 and in Figure 8 (center). Testing with the NoisyTestDataSET while the MLP was trained with CleanDataSET shows (as expected) a large increase in the average error.
More robust results with respect to noise can be obtained by training again the best MLP on the full NoisyTrainingDataSET. We do not perform an additional hyperparameter tuning, but simply repeat the training stage on the same MLP. In this case, the larger dataset and the presence of noise induce an error when tested with the NoisyTestDataSET of around 5.92%. This error can, again, be improved by removing the average as we did before, leading to a final error which is very close to the clean case (4.96% vs 4.67%). We show three doping reconstructions for this case in the 25, 50, and 75 percentile in Figure 9.
In summary, introducing a nonlinear MLP model has reduced the error roughly by a factor of two compared to the linear least square model. The average MLP is just below 5% for clean test data or noisy but properly average adjusted test data.

Residual neural network (ResNet).
In 2015, Kaiming et al. [13] introduced ResNets, feedforward convolutional neural network models, which, since then, have been widely used to solve different kinds of problems. One of the biggest advantages of ResNets (or, generally, of convolutional neural networks) over the multilayer perceptrons is their ability to learn 2D images. Interpreting the doping profiles in 2D probing region as an image, this is precisely what we need to achieve, too. Moreover, the proportionality shown in Equation (4.6), valid under additional assumptions, suggests that the photovoltage signal and the doping profile may be related by convolutions. Even assuming now a nonlinear relationship described via ResNets, we may still wish to exploit a convolutive structure to encompass the relationship in Equation (4.6). Our reference implementation for residual neural networks (ResNet) is described in [13], with some important differences regarding the structure of the network.
First, in our LPS setup from Subsection 4.2, we consider one-dimensional data vectors C and u instead of two-dimensional images. The dimension is particularly relevant to choose the size of the network: In ResNets, the downsampling blocks reduce by a factor of two the size of the input, by first halving the size of the input in each direction (reducing the spatial indices by a factor four) and then doubling the number of channels. In our setup, this dimension reduction associated to the downsampling block does not happen because in the first step we only have one dimension to divide by two, and we end up with the same number of neurons after doubling the numbers of channels. Moreover, while the ImageNet database used to train the network described in [13] contains about one million images, our datasets are significantly smaller. We expect to face some overfitting and therefore we aim for a model with fewer parameters than the one described in [13]. Finally, we solve a regression problem instead of a classification one. The major consequences for our model are that we cannot use drop-out layers to reduce overfitting, and it is unclear if there is any benefit in using batch normalization layers. We keep the batch normalization layers because we expect the statistical distribution of our batches to be similar to each other, and the statistical distribution of the error for a simpler network (the MLP) shows an improvement in the performance of the model with noise when removing the average, suggesting that batch normalization layers are at least not harmful. We build several ResNet models using PyTorch ( [32]), and develop a strategy to find the best ones by adapting the model described in [13], and using the insight obtained from the MLP model presented in Subsection 4.5. Our best multilayer perceptron had 160,950 parameters that had to be optimized. We try to keep the number of parameters of our ResNet around the same order of magnitude. Of course, this means that we need to drastically simplify the model developed in [13], to find a model whose number of parameters is compatible with the size of our

CleanTrainingDataSET.
Since training a ResNet is significantly more expensive compared to training a multilayer perceptron, we use a two-step strategy to reduce the computational cost of the hyperparameter tuning. The possible configurations of the ResNet that we allow in our models are detailed in Appendix B, and correspond to a total of 48 configurations for the gate, 9 configurations for the encoder, and 18 configurations for the decoder, for a total of 7776 possible configurations.  In the first step of our hyperparameter tuning, we restrict the parameters for the SGD optimizer to a fixed float number chosen in the interval [5 × 10 −3 , 1 × 10 −1 ] for the learning rate, and choose a batch size in the set {64, 128, 256, 384, 512}. We fix also the gradient clipping to 1 and the weight decay to 0. During this step, we randomly sampled 300 ResNet configurations among the 7776 generated according to Appendix B, and for each one of them, we perform six trainings using different set of parameters for the SGD optimizer. Again, we rely on the ASHA algorithm as we did in Subsection 4.5. We compare the performances of the models after 5000 epochs of the SGD algorithm, and retain only four of the candidate ResNet models selected by the ASHA algorithm, which are described in Table 1  In the second stage of the hyperparameter tuning, we keep the four structures of the ResNet models fixed, and enlarge the search space in the optimization parameters by defining a new set of admissible parameters: we choose the learning rate in the interval [10 −3 , 10 −1 ], the batch size in the set {32, 64, 128, 256, 384, 512}, the gradient clipping in {10 −2 , 10 −1 , 1}, and we introduce some weight decay, with decay parameters chosen in the set {10 −4 , 10 −3 , 10 −2 , 10 −1 , 1}.
For each one of the four models, we perform 200 different trainings applying different parameters of the SGD optimizer, and finally select a winner, corresponding to RN3 in Table 1.
In Figure 11 we present the statistical distribution of the error for ResNet RN3 selected with the strategy outlined above. We observe that the average error we obtain using the CleanDataSET for both training and testing is around 5.47% (see Table 2 for the details). This number is slightly worse compared to the MLP case introduced in Subsection 4.5. In particular, we observe that the ResNet RN3 trained with the CleanDataSET is more sensitive to noise (see Figure 11, center figure) when compared to the MLP (29.7% vs 21.2%). This difference decreases when removing the average, but for the ResNet remains surprisingly worse than the least square model (16% vs 15.5%).
When we keep the structure of RN3, and train the network on the NoisyTrainingDataSET, we obtain the results shown in Figure 11 on the right. An example of reconstruction obtained with RN3 with a sample from the NoisyTestDataSET is shown in Figure 12. In this case the average error on the NoisyTestDataSET is around 6.72%, which is slightly worse than the result obtained with the MLP network.

Summary and outlook.
The aim of this paper was to properly model and numerically solve general ill-posed inverse photovoltage technologies where measured photovoltage signals are used to reconstruct local doping fluctuations in a semiconductor crystal. The underlying charge transport model is based on the van Roosbroeck system as well as Ohm's law. We presented three different data-driven approaches to solve a physically relevant local inverse problem, namely via least squares, multilayer perceptrons, and residual neural networks. The methods were trained on synthetic datasets (pairs of discrete doping profiles and corresponding photovoltage signals at different illumination positions) which are generated by efficient physics-preserving finite volume solutions of the forward problem. While the linear least square method yields an average absolute ∞ error of 9.3%, the nonlinear networks roughly halve this error to 4.7% and 5.5%, respectively, after optimizing relevant hyperparameters. Our method turned out to be robust with respect to noise, provided that training is repeated with larger, noisy, datasets. In this case, the error is around 9.8%, 5%, and 6.7% respectively. Removing the average doping value from the data was more important to reduce the testing error for clean datasets than for noisy ones. The datasets, python codes, and resulting trained networks are available in the repository [34].
Future research may go into different directions: the numerical simulation served as a proof of principle and was limited to a 2D version of the photovoltage problem. For the 3D problem, efficient data generation strategies need to be designed. Furthermore, it is clear that different doping profiles may correspond to similar signals. On the one hand, this is intrinsically dependent on the technology that we apply to recover the signal: we hit the crystal with a laser that has a   finite laser spot radius (in our case, around 20µm), and we cannot expect to resolve oscillations with smaller wave lengths. On the other hand, our models do not have a way to distinguish between two doping profiles that deliver the same signal, but are exposed, during training, to many such cases. In Figure 13, for example, we show two examples of doping profiles that have errors on the far right regions of the histograms in Figure 8 and Figure 11. The figure shows how in some cases, the inverse problem is oblivious to higher oscillations, and both neural network types will have cases in which they return an answer which is either too smooth w.r.t. the expected result, (left panel of Figure 13) or too oscillatory (right panel of Figure 13). A possible focus of future research could be to study how to resolve this ambiguity, which is intrinsic to the ill-posedness of the discrete local photovoltage inverse problem. Finally, we aim to extend our data-driven approach to opto-accoustic imaging, replacing the charge transport model with appropriate thermal expansion models. factor two, while the fixed channel do not. In this case the shortcut operation of both blocks is a convolutional layer with kernel size equal to one and stride equal to two (basic block), or kernel size equal to two and stride equal to two (fixed channel block), followed by a normalization block. The convolutional layer of the shortcut of the fixed channel block forbids communication between channels, i.e., each element of the output tensor depends only on the values of the elements of the input tensor that share the same index for the channel (or, in other words, using PyTorch we impose that the number of groups of the convolutional layer is equal to the number of channels). For what concerns the computation of the output (and not the shortcut), the dimensional reduction is obtained by setting the stride of the first convolutional layer to 2.
When downsampling is not active, the shortcut operation of both blocks is the identity tensor. Therefore, the only free parameters that we have left for our encoder are the number of blocks, their kind, and a downsampling flag for each block. Our encoders consist of one, two, or three blocks of the same kind. For encoders made of basic blocks, we allow two different configurations of the downsample flag: true for all blocks or false for all blocks. For fixed channel blocks, we always set the downsample flag to true. We have a total of 6 configurations for the encoders with the basic blocks and 3 configurations that use the fixed channels blocks, for a total of 9 possible configurations.
Decoder. This is essentially a multilayer perceptron, and we allow 1 or 2 hidden layers. If we choose a configuration with 1 layer, the size of this layer could be 100, 150 or 200. With 2 layers, there is a total of 15 possible configurations obtained by choosing the size of the first layer in {100, 150, 200, 250, 300} and the second one in {100, 150, 200}. In total, we have therefore 18 different configurations for the decoder.