Mathematical Modeling of the Limiting Current Density from Diffusion-Reaction Systems

The limiting current density is one of to the most important indicators in electroplating for the maximal current density from which a metal can be deposited effectively from an electrolyte. Hence, it is an indicator of the maximal deposition speed and the homogeneity of the thickness of the deposited metal layer. For these reasons, a major interest in the limiting current density is given in practical applications. Usually, the limiting current density is determined via measurements. In this article, a simple model to compute the limiting current density is presented, basing on a system of diffusion–reaction equations in one spatial dimension. Although the model formulations need many assumptions, it is of special interest for screenings, as well as for comparative work, and could easily be adjusted to measurements.


Introduction
In electrochemisty and chemical engineering, the limiting current density indicates the state at which the maximum metal deposition can take place in a galvanic deposition, see [1,2]. For this reason in practical applications, the limiting current density is often used as an indicator for the maximal deposition speed during electroplating. Additionally, note that due to Faraday's laws, see [1][2][3][4], the deposited mass of metal is proportional to the applied current on the cathode surface, which differs on complex geometries, see [5].
With the limiting current density as an indicator of the upper bound for the localized metal deposition speed, it is an indicator for the attainable homogeneity of the deposited metal layer as well.
The limiting current density of an electrolyte is most commonly determined with electrochemical measurements, such as current-potential curves, see, for example, [1,3]. With these means, electrolyte design via the limiting current density can be very time consuming and expensive in terms of needed material, especially if screenings are taken into account. For this reason, a model is needed, which can be used as an additional tool to support the screening and hence needs to be relatively efficient to evaluate.
In the literature, both globalized models, see [1,3,5,6], considering the processes in a galvanic bath with descriptions of static electric fields at their core, and localized models, considering the processes directly at the electrodes concentrating on diffusion-reaction formulations, see [2,3,[7][8][9][10][11], exist and could give hints for the limiting current density. Stateof-the-art approaches to model the limiting current density, as discussed in [12], use the assumption of thermodynamic equilibrium in the bath and use the corresponding species distributions. However, those approaches omit the processes in the diffusion boundary layer. Hence in the following work, a model to compute the limiting current density by combining the idea description of the species concentrations by diffusion reaction equations with the idea from [12] and a model adaption from the model in [13], assuming thermodynamic equilibrium in the bath.
In Section 2, the mathematical model for the limiting current density will be derived, basing on a three-component system of diffusion-reaction equations describing the species concentrations in the diffusion boundary layer and will be treated with the numerical methodology described in [14]. A method to evaluate the thermodynamic equilibria is stated in [15].
In Section 3, the model will be applied to a Cu-CO 3 electrolyte and demonstrates the applicability of the methodology, which is described in [14,15], to a real-life problem and will be compared to a well-known model with no reactions in the diffusion boundary layer, described in [13].
Lastly, in Section 4, the results of the article will be summarized and discussed. Additionally, an outlook for future works will be given.

Mathematical Modeling of the Limiting Current Density
This section is devoted to the mathematical modeling of the limiting current density. This model strongly relies on the description of the localized species concentration distributions in the diffusion boundary layer at the cathode. The description starts with the general setup of the physical experiments.
Assume that there are a cathode and an anode positioned in an aqueous electrolyte consisting of a metal M, a ligand L and complexes K κ , for 1 ≤ κ ≤ R ∈ N, produced from the following reactions where l κ is the stoichiometric index of the ligand L, m κ is the stoichiometric index of the metal M, h κ is the stoichiometric index of the protons H + and o κ is the stoichiometric index of the hydroxide-ions OH − . Additionally, assume that in the bath, except for the diffusion boundary layers at the electrodes, the species-concentrations c S,eq , for S ∈ {M, L, K 1 , . . . , K R } are in thermodynamic equilibrium. For further discussion, note that numbers and functions are commonly indexed with corresponding species. To shorten the notation, where it is useful, in this article species are replaced by the index S and a list of species that can be inserted as S.
As discussed in [15,16], the concentrations in thermodynamic equilibrium fulfill the following mass conservation laws where m S , for S ∈ {L, M, OH, H}, denotes the total masses of the species S in the bath.
To complete the mathematical system describing the thermodynamic equilibrium, as in [15], the definition of the complex-formation constant 0 < β κ of the κ-th reaction in (1), see [16], is used: As discussed in [15], the Equations (2a)-(2e) describe a model for the concentrations in thermodynamic equilibrium in the bath and are uniquely solvable for non-doubling definitions of complexes.
To describe the cathodic processes it is assumed, as in [7,8], that, locally, at the cathode in the diffusion boundary layer of the width δ X , the species distributions are not determined by thermodynamic equilibrium but by diffusion and the following single reaction: Additionally, assume that the only geometrical property due to which the concentration c S , with S ∈ {M, L, K}, changes is the distance to the cathode. Without loss of generality, assume that the cathode surface is positioned at x = 0, hence the domain in which the discussion takes place is Ω =]0, δ X [. As defined in [1], the limiting current density is the maximal current density at which the deposition of the metal M can take place. Hence when applying the limiting current density, maximal consumption of the metal M takes place. Thus, the metal-concentration at the cathode vanishes in the kinetic equilibrium. Hence, it suffices to consider static diffusion-reaction laws to compute the limiting current density and set c M (0) = 0 as boundary condition. Additionally, suppose that no species transport of the species L or K takes place over the cathode surface, i.e., d dx c S (0) = 0, for S ∈ {L, K}. With the discussion and assumptions above the concentrations c L , c M and c K the species M, L resp. K can be described by the following boundary-value problem, see [2]: Note that the problem formulation (4a)-(4e) does not include the positivity of the species concentration, which will lead to negative species distributions, as discussed in [14]. A common strategy to include the positivity of the concentrations into a model formulation is the usage of restrained minimization problems, see [17]. A reformulation of (4a)-(4e) was undertaken in [14] and relies on the introduction of new variables p S = d dx u S , for all S ∈ {M, L, K}. The introduction of the new variables yields the following mathematical problem.
Find u S , p S S∈{L,M,K} ∈ H 1 (Ω) ×6 such that the following equations hold true: Note that with H 1 (Ω) the vector space of once weakly differentiable L 2 functions is denoted, cf. [18,19], i.e., the functions for which partial integration is well defined.
Assuming the mathematical problem (5a)-(5f) is uniquely solvable it is equivalent to the following least-squares problem: Find x ∈ X = H 1 (Ω) ×6 such that the following equation holds true: where the nonlinear least-squares functional LS ∶ X → R for all y = (v S , q S ) S∈{M,L,K} ∈ X is given by and is obtained by taking the squares of the L 2 (Ω) norm of the Equations (5a)-(5d) and taking the absolute norms of the Equations (5e)-(5f). Summing the terms, one obtains the corresponding least-squares functional and sees that a unique minimizer, (6), is also a solution of (5a)-(5f).
The solution of (6) allows, in general, negative species concentrations 0 < c S . This makes a model adaption necessary; the model adaption is the reformulation of (6) into an obstacle problem as done in [17]. This is generally done by restraining the set of functions, over which one wishes to minimize.
The resulting problem formulation reads as follows: Find x ∈ X = H 1 (Ω) ×6 such that the following equation holds true: From Faraday's law, see [1,2]; it follows that the limiting current density j lim is given by: where F is the Faraday's constant and z is the valence of the metal M and (c S , p S ) ∈ X is the solution to (7). With the representation (8) of the limiting current density j lim . At the end of this section some remarks must be made.

Remark 1.
Note the following remarks: (i) First note that the systems of ordinary differential Equations (4c)-(4e) and (5c)-(5f) is due to the theorem of Picard-Lindelölf-see [20]-which is uniquely solvable. (ii) A robust and efficient methodology to compute approximations of (7) is described in [14], which is based on a combination of augmented Lagrangian methods, conforming adaptive FEM, quasi-Newtonian methods and homotopy methods. Although the methodology is quite expensive in terms of CPU time, it is efficient in comparison to other state-of-the-art methodologies, cf. [14]. (iii) For known RHS m L , m M , m H , and m OH , the model Equations (2c)-(2e) are uniquely solvable and a stable and robust methodology is described in [15]. (iv) The model equations for the species distribution in the diffusion boundary layer (4a)-(4e) can easily be extended, in theory, to more reactions with analogous boundary conditions, with concentrations in thermodynamic equilibrium at the bath-end of the boundary layer and no material transport over the cathode surface for the species except the free metal concentration.
Additionally the definition (8) stays basically the same. The theoretical work in this article reduces to the one major reaction with the biggest effect, since, as known to the authors, no effective algorithm for the numerical treatment of the positivity of the respective species is given. (v) Note that the thickness δ X , which is an estimated value in many cases, has a strong influence on the calculated value j lim , due to its influence on the solution of (7). (vi) Assuming that no reaction takes place in the diffusion boundary layer and the only process which takes place is diffusion of the metal species M. Assuming the same boundary conditions for c M as in (4d) for the diffusion law directly implies, with the unique solvabillity of (9) and the assumed boundary conditions, the following representation of c M ∶]0, δ X [→ R: In this case, the limiting current density reduces to This formula coincides with the formula given in [13]. (vii) The exact representation (11) yields a first formula for the calibration of δ X , which is needed, see (iv). Assuming that a measured limiting current density j lim,M is given, then δ X can be approximated via

Application to the Cu Deposition from a Cu-CO 3 Electrolyte
In this section, some practical examples to the computation of the limiting current density will be discussed. The necessary parameters, i.e., reactions, diffusion coefficients, rate coefficients and stability constants, for the Cu-CO 3 electrolyte discussed in this section can be found in [10]. Additionally, a discussion of the limiting current density over the pH value is discussed.
In Table 1, the species, reactions and, complex formation constants to describe the thermodynamic equilibrium in the bath, see Equations (2a)-(2e), are given. Additionally, the diffusion coefficients and rate constants for the corresponding species and reactions in aqueous media are given. In the following, assume that, in the diffusion boundary layer, only the reaction takes place. In [10] the thickness of the diffusion boundary layer 0 < δ X is considered to have the length δ X = 21 µm and a valence of z = 2 for the Cu ions is assumed. The methodology used for the simulation of the corresponding diffusion-reaction laws (see problem (7)) is described in [14]. The numerical methodology described in [14] includes a combination of augmented Lagrangian methods, a conforming lowest order discretization, adaptive mesh refinement, homotopy methods and quasi-Newtonian methods. The estimator is given by the evaluation of the corresponding least-squares functional in the discrete approximation x h ∈ X h . In this section, four sets of numerical experiments will be discussed. The numerical examples are defined by fixing m Cu as m Cu = 0.5 mol l for all numerical experiments, and a definition a simple variation of the total masses m CO 3 = 0.5 mol l , m CO 3 = 0.75 mol l and a variation of the thickness of the diffusion boundary layer δ X = 21 µm, δ X = 31 µm. Additionally, the total masses m H and m OH were chosen such that a pH interval [2,12] is covered for the discussion of the limiting current density.
In Figure 1, the limiting current density distributions over the pH value for δ X = 21 µm and m CO 3 = 0.5 mol l (a) and for δ X = 21 µm and m CO 3 = 0.75 mol l (b) are shown. The limiting current density for δ X = 31 µm and m CO 3 = 0.5 mol l (c) and for δ X = 31 µm and m CO 3 = 0.75 mol l (d) are shown as well. As it can be seen in the figure, the current-density distributions for m CO 3 = 0.75 mol l are lower than for m CO 3 = 0.5 mol l . This can be expected, since the more CO 3 is in the system the higher is the mass of metal ions bound to the CuCO 3 complex and, hence, the free metal concentration c Cu is lower. A second observation is that the current-density distribution is strongly dependent on the assumed thickness δ X of the diffusion boundary layer. The last observation is that the limiting current density seems to decrease as the pH value increases. This observation corresponds with the expectations since the free metal concentration in thermodynamic equilibrium decreases as the pH value increases.  In Figure 2, the species concentrations for the Cu species, CO 3 species and the CuCO 3 species for m CO 3 = 0.5 mol l and δ X = 21 µm at a pH value of 2 are shown. As can be seen in the figure, the Cu concentration decreases with the distance to the cathode. The same behavior can be seen for the complex concentration c CuCO 3 . This behavior can immediately be explained by the reaction kinetics, i.e., the mass of free Cu ions decreases by the deposition at the cathode. This effects the mass of free Cu ions that can be bound in the complex as well. As a consequence, the ligand concentration increases as the distance to the cathode decreases. In Figure 3, the Cu concentrations c Cu are depicted at a pH value of 2. In subfigure (a), the concentration profiles c Cu of the Cu species for m CO 3 = 0.5 mol l , 0.75 mol l and δ X = 21 µm are shown. In subfigure (b), the concentration profiles c Cu of the Cu-specie for m CO 3 = 0.5 mol l , 0.75 mol l and δ X = 31 µm are shown. As seen in both subfigures, the concentrations decrease as m CO 3 increases. In this setting, this means that the more CO 3 is in the system the more Cu is bound in the complex CuCO 3 , which makes this result expected.  In Figure 4, the estimated error of the approximate solutions x h are shown. In this context, the error was estimated by the evaluation of the least-squares functional LS(x h ).
As it can be seen in Figure 4, the evaluation of the error estimator experimental setups and pH values. As can be seen directly, the simulations used for the experimental evaluations are in the convergent part of the graphs and have, in comparison to the initial errors, sufficient low evaluations of the error estimator such that the simulations can be considered reliable. In Table 2 a comparison of evaluations of the limiting current density (LCD) w.r.t. the model described in [13] and derived in this paper is undertaken. The evaluations vary over pH values, the discussed mass differences of CO 3 m CO 3 = 0.5 mol l and m CO 3 = 0.75 mol l and the thickness δ X of the Nernst diffusion layer δ X = 21 µm and δ X = 31 µm. As seen in the table the following observations can be made. Both models predict, in each case, a decrease of the limiting current density as the pH value increases. Furthermore it can be seen that, in both models, the limiting current densities for δ X = 21 µm tend to be higher than for δ X = 31 µm. Additionally it can be seen that, as expected, the limiting current density is in both models higher, when considering higher ligand masses. Table 2. List of computed limiting current densities (LCD) on the basis of the model described in [13] and the manuscript described in Section 2. Comparing the values for the new model and the model described in [13] in Table 2, one sees that the computed values of the LCD from the new model are a factor of 10 3 higher than the LCD from the old model. Exceptions from this rule are the values for δ X ) = 31 µm and m CO 3 = 0.75 mol l , where the limiting current density from the old model is at the factor 10 3 over the new model. To explain this asymmetric behavior, note that the formula for the LCD in the old model (see [13] and (11)) is derived from an exact and linear solution to the second Fickian law representing the species distribution c M , which implies a linear dependency of the LCD of 1 δ X and c M,eq . In contrast to that the species distribution, c M is highly nonlinear and, thus, the distribution of the limiting current density over the pH value can be expected to be nonlinear. This can and will explain the asymmetry addressed.

pH-Value
While both models are dependent on the pH value, the total mass of the ligand m CO 3 and the thickness of the boundary layer δ X , the variance of the values with the new model seems to be higher and to underly a more complex nonlinear law than described in [13]. Additionally, it can be seen that, in most cases, the predicted limiting current density from the new model is higher than the respective value from the model described in [13].
In total, the new model behaves as expected. Although a strong dependency of the limiting current density of the thickness δ X of the diffusion boundary layer is shown, the experiments also show that the model can be used to study the pH dependencies of the limiting current density of an electrolyte and the recipe of the electrolyte can be studied as well. Furthermore, the model can be applied to determine possible alternatives for toxic ligands.

Discussion and Conclusions
In this paper, namely in Section 2, an abstract model for the limiting current density was derived. The model described in Section 2 assumes thermodynamic equilibrium of a broad system of reactions in the bath (1), described by the Equations (2a)-(2e). Additionally, it is assumed that, in the diffusion boundary layer at the cathode, only one reaction, namely (3), takes place. Additionally it is assumed that the species concentrations w.r.t. the metal M, the ligand L and the complex K obey the diffusion in the boundary layer. The corresponding physical behavior was described by a stationary system of diffusion and reaction-diffusion equations. Through a convenient definition of the boundary conditions and a reformulation of the arising problem into an obstacle problem to treat possible negative species concentrations, it was possible to compute the limiting current density j lim from the metal concentration c M through (8).
In Section 3, an electrolyte was numerically studied. As a result of the numerical experiments, it can be seen directly that the limiting current density is strongly dependent on the mass input m L and m M as well as the assumed thickness δ X of the corresponding boundary layer.
Although the model formulation for the limiting current density is relatively simple, especially only using the language of ODEs (ordinary differential equations), it is strongly biased. For instance, the system of reactions in the bath, which is assumed to be in thermodynamic equilibrium, may and will be much larger than the single reaction assumed to take place in the diffusion boundary layer, but can be expected to be more precise than the model described in [13], since first the thermodynamic equilibrium at the boundary of the Nernst diffusion layer is introduced and a reaction is introduced in the layer domain; hence, missing reactions are respected.
As discussed in Section 3, the computed value of the limiting current density strongly depends on the assumed thickness of the diffusion boundary layer δ X . Similar results, concerning the dependency on δ X , are presented in [8]. Besides the formula given in Remark 1 there are approaches to computing the thickness, as in [10], relying on the currents in the layer but having its own biases. A convenient way to determine the thickness of the diffusion layer would be a direct measurement of the boundary layer. A method to measure it is described in [21], which relies on the measurement of the changes in pH value under the change in time. Further measurements for laminar and turbulent flows of electrolyte are described in [22].
Another bias, which it is necessary to discuss, is the assumption of the thermodynamic equilibrium in the bath. As discussed in [1,4,5,23], the transport processes in the bulk electrolyte are dominated by diffusion, migration, convection and reactions. Assuming, as in [7,8,11,12], that the concentration gradients are low, a concentration gradient will normally exist and, hence, thermodynamic equilibrium will not be fulfilled. Additionally, as it is commonly known, see [24][25][26], the transport processes are strongly geometrydependent and hence, gradients of the species distributions are also geometry-dependent. In particular, it is typically assumed that the diffusion boundary layer lies within a laminar boundary layer, the thickness of which will differ along the corresponding cathode surface, which will bound the thickness of the diffusion boundary layer due to the geometry dependence of the mechanical flow of the electrolyte (see [26]).
Beyond the biases discussed above, the model has various advantages. For instance the model is relatively easy, and it can be evaluated with a relatively small numerical cost. Even though deviations between the experiment and the model prediction are expected, the model can be used for a more general, qualitative comparison. For example, as discussed in Section 3, electrolytes can be studied within a variation of the pH value. Additionally, the software can be used to compare electrolyte mixtures wherein the ratio of m L to m M is varied. Additionally, the model can be used to determine possible substitutes for toxic ligands. The last point is of considerable importance in designing, for example, REACHcompliant electrolytes. In those cases, the model will be sufficient.
In summary, a model was introduced, which, on the one hand, is highly biased, and has, on the other great value in the comparative development of electrolytes. For the formulation of a quantitative, more precise model, more research is needed. The simplest approach to defining such a model is to extend the model to a geometrically and specifically known setup, such as the description of the complete cell, only considering the distance from anode and cathode as the important geometrical feature. This geometrical adjustment needs the inclusion of migration and further reactions assumed to be not in thermodynamic equilibrium. Funding: The publication of this article was funded by Chemnitz University of Technology.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: subject to w.r.t. with respect to resp. respective LCD limiting current density ODE ordinary differential equation