General Approach to Model the Surface Charge Induced by Multiple Surface Chemical Reactions in Potentiometric FET Sensors

We propose a general methodology to calculate the individual sensitivity and the cross-sensitivities of potentiometric sensor devices (e.g., ion sensitive FETs (ISFETs), CHEMFETs) with an arbitrary number of non-interacting receptors binding to ionic species or analytes in the electrolyte. The surface charge generated at the (bare or functionalized) interface with the electrolyte is described by the Poisson equation coupled to a linear system of equations for each type of receptor, where the unknowns are the fractions of sites binding with a given ion/analyte. Our general model encompasses in a unique framework a few simple special cases so far separately reported in the literature and provides for them closed-form expressions of the average site occupation probability. Detailed procedural description of the usage and benefits of the model is shown for specific cases with concurring surface chemical reactions.

because equally charged ions or analytes cannot be distinguished from each other. Achieving high selectivity is, in fact, a challenging goal when developing ion-sensitive devices. Detailed understanding of the FET transduction mechanisms is important to achieve a quantitative description and prediction of the sensor response, selectivity, and cross-sensitivity to multiple charged analytes useful to interpret experiments and develop an optimized device and circuit designs [13]- [15]. Ad hoc models have been successfully used to describe the experimental observations [16]. However, a general model has not been developed so far amenable to incorporate all special cases in a comprehensive theory. The approach we propose in the following sections generalizes these ad hoc models to calculate the surface charge at solid/electrolyte interfaces. Simple graph charts of the surface reactions yield in a straightforward manner the equations to compute the expected net surface charge. The theory behind the model is given in Section II, where case studies are also illustrated. Section III develops a 1-D simulation framework to demonstrate the main features of our approach. Finally, in Section IV, we validate the model by reproducing a few sensitivity calculations and experiments from the literature. This article expands significantly over our conference presentation in [17] by reporting the full details of the theory, new closed form equations, analysis of more relevant cases, and demonstration of a calibration procedure.

A. Starting From a Simple Example: The Langmuir Reaction
To introduce the method and notation we start with a uniform and perfectly flat interface between a solid and a solution and a uniform layer of ligands. Let f o and f f = 1 − f o be the average fraction of ligands that are occupied by (respectively, free of) a given analyte A having volume concentration [ A]. In this case, adsorption and desorption mechanisms are generally addressed with the simple Langmuir adsorption model [18]. The adsorption rate R a (similar to trapping/detrapping mechanisms in semiconductor physics [19]) is assumed proportional to the concentration of analytes in the solution, their thermal velocity v th and the capture cross section ; that is, R a = [A]v th . Conversely, the desorption rate R d depends on the binding lifetime τ A as R d = 1/τ A . The time-dependent occupation probability, f o (t), is determined balancing the net flux of analytes to the surface Such a system can be represented by the transition diagram in Fig. 1(a). In the steady state, the time average occupancy of each state is easily computed setting d f o /dt = 0 and rearranging (1) to obtain According to (2), we recover a Langmuir-type reaction [18] where the only parameters are the analyte concentration and the dissociation constant. Equivalently, we can simplify the transition diagram into a graph chart, where the coefficient of the branch between two nodes (states) is the analyte concentration divided by the dissociation constant, see Fig. 1(b). In fact, the mathematical relationship expressed by the graph chart (2). It is thus possible to derive the steadystate equilibrium solution from the graph chart coupled to the normalization condition.

B. General Formulation
Let us consider a generic set of binding reactions between surface sites/ligands and electrolyte ions/analytes. We assume that the site or ligand can bind in sequence with an arbitrary number of ions/analytes of different species (see Fig. 2). Our analysis starts by keeping track of all the possible complexation forms, assigning to each a state and its net signed number of elementary charges, z i . For a given site, we can then identify the states and the reactions that transform one state into another one. In fact, each reaction corresponds to the binding/unbinding of a given ion/analyte from one complexation state to another. As a result, N −1 reactions will create links between N states (as the empty state, i.e., ligand not bound to any ion, is also considered). The graph representation of a generic reaction will be composed of N − 1 arrows (each representing a specific reaction) connecting the N nodes. Similar to Fig. 1(b), each arrow has a coefficient given by the ratio between the concentration of the ion/analyte and the corresponding dissociation constant. For example, if the ligand has probability f i to be in state "i," the probability to find it linked to species A, hence in the new state " j," is where K i j is the dissociation constant of the adsorption reaction. This leads to N−1 equations expressing the dependencies between the state probabilities of two consecutive states. Note that for the sake of clarity and of a simple notation, all the chemical reaction constants used throughout this manuscript are dissociation constants, consistently with the majority of the experimental data shown in Section IV. As for the previous example of the Langmuir reaction, the missing closure equation is given by the normalization condition N 1 f i = 1. Hence, a set of N equations and N unknowns is obtained which can be easily written in a matrix form M · f = (0 · · · · · · · 0 1) T (4) where f = f 1 · · · · · · · f N T and the N × N matrix M takes the form The firsts N − 1 rows represent the equilibrium state of the chemical reactions while the last row (a vector of ones) is the normalization condition. A generic reaction is highlighted in red: By construction, if the n th reaction defines the transition of the ligand between state i and state j , then the i th column of the n th row contains the transition coefficient according to the graph while the j th column contains the constant −1. The matrix expresses the relationships among the unknown state probabilities, which can be computed as where b = 0 · · · · · · · 0 1 T . Once the state probabilities f i are known, it is straightforward to compute the surface charge. In fact, by assuming uniformly distributed and independent binding sites we can calculate the fraction of sites that, at equilibrium, are in state i and then the net surface charge per unit area Q S at the sensing layer at equilibrium where q is the absolute value of the elementary charge, N S is the total number of sites per unit area and z i is the net signed number of elementary charges (which does not need to be integer) of state i . The Q S should be properly taken into account when solving the Poisson-Boltzmann equations for the electrostatics as explained in Section III.

C. Two Relevant Special Cases
A direct generalization of the Langmuir process is given by the chained reaction, where each site undergoes consecutive reactions in which chemical species are sequentially adsorbed.
A second relevant case occurs when one binding process excludes further adsorptions, but the same site can accept different ions/analytes which then compete for binding to the same site. These two important cases can be addressed analytically and the graph properties lead to compute each state probability as outlined below.
1) Chained Reactions: Let us consider the sequence of binding reactions, illustrated in the simplified graph chart of Fig. 3. The mathematical representation of this graph leads to a set of N − 1 equations, each of them linking the state probability f i+1 associated with a state with the one associated with the previous state f i , multiplied by the respective branch coefficient, i.e., After few algebraic steps, the expression of the fractional occurrence f i is given as 2) Interference by Multiple Ions/Analytes: A second relevant case that can be addressed analytically is the multiple interference of N −1 different ions/analytes that compete to bind with a site/ligand L. The graph chart is given by a branched diagram (see Fig. 4). As for the previous case, we can derive the N × N matrix from the graph including the normalization of the f i functions. By denoting with suffix "0" the state corresponding to the unoccupied ligand, we can write ⎛ Simple algebra leads us to compute the state probabilities f i as For example, consider two interfering ions A and B binding to the same neutral ligand L. If the respective chemical reactions have dissociation constants K A and K B and lead to a net signed number of elementary charges z A and z B , we obtain the surface charge density Note that (12) is the same derived in [20] for the competition between proteins and ions in linking to a recognition molecule.

D. More General Cases
The method allows users to describe much more complex sets of reactions so far as one simple rule is satisfied: In constructing the graph chart, two arbitrary states must be connected by only one path. Consider the example in Fig. 2. From left to right, the ligand L binds to ion C and then ion A aggregates to the complex leading to the state LC A. In this case, the process can also happen stepping from L to L A and then L AC and the complex L AC has no chemical distinction with respect to LC A. This would lead the graph chart to show two different paths going from L to LC A = L AC, giving more equations than unknowns. However, in modeling such situations, one should construct the graph so that separate paths stemming from the same state always lead to different states even if, chemically speaking, these final states are equivalent (see Fig. 2 on the right). By considering LC A and L AC as two separate but indistinguishable states, the respective state probabilities f LC A and f L AC can be derived. The indistinguishable states LC A and L AC will then both contribute to the interface charge. Fig. 5 provides another example of chemical reaction other than the linking of a single ion/analyte to a site. In fact, the graph describes the general reaction L A − + B + L B + + A − where the factor 1/K A K B represents the dissociation constant of the reaction and the ion concentrations in the branch expressions are at the numerator or at the denominator depending on whether the ions are absorbed or released by the reaction.

E. Consistency With Known Models
Our methodology embraces in a general framework several models proposed in the literature to compute the surface charge at equilibrium. Here, we present two cases relevant for ISFETs.
1) Site-Binding (SB) Model: The SB model [10] is widely used to describe potentiometric pH sensors. According to this model the hydroxyl groups in the outermost atomic layers of several metal oxides behave as amphoteric sites that interact with hydrogen ions in the solution. A combination of protonation and deprotonation reactions takes place and leads to the creation of a net surface charge. By denoting the generic metal with M and its respective hydroxyl state as MOH, we identify the protonated state as MOH + 2 and the deprotonated one as MO − . Alternatively, these states can be seen as the sequential adsorption of protons, starting from MO − and considering two protonation reactions. As a result, we identify three states and two chemical reactions, which entails a simple chained-like graph as illustrated in Fig. 6. The Fig. 6 is the proton concentration at the surface, K a and K b are the commonly adopted acid and basic dissociation constants [10] while f i and z i have their usual meaning. In the matrix form, we obtain ⎛ Graph chart of the SB model [10]. Two protonation reactions (arrows) transform the negative MO − state to the positive MOH + ¾ state. Each protonation reaction involves the binding of a proton and is ruled by the dissociation constant K a or K b . Fig. 7. Graph chart for the mSB model [11]. Compared to Fig. 6, the chained reaction involves a new state labeled MOHCl − which describes the adsorption of chloride ions and the expulsion of a proton.
By solving for f we can estimate the surface charge at equilibrium as where N S is the number of sites per unit area. The expression in (14) is identical to the one of the original work in [10] and to (9) if one notes that the division of all terms by K a K b results in the sum of f 1 and f 3 .
2) Modified SB (mSB) Model: The modified version of the SB model proposed in [11] to account for the influence of chloride ions can also be derived using our approach. In this case, an additional reaction involving the adsorption of chloride ions and expulsion of a proton MOH + 2 + Cl − → MOHCl − + H + is added to the previous case. Using the notation of [11] for the reaction constants and the rules given in Section II-D the result (see Fig. 7) is characterized by four states (MO − , MOH, MOH + 2 , and MOHCl − ), two protonation reactions and a third new reaction describing the surface complexation of chloride ions. Translating into the matrix form, we obtain ⎛ Once the states' probabilities are known, we can calculate the surface charge (16), as shown at the bottom of the next page, in agreement with the original work [11]. Note that, following [11], we assume that the reaction with the Cl − anion involves emission of a proton to the bulk of the electrolyte solution; hence the concentration is denoted as [H + B ], differently from the [H + S ] involved in the protonation of the MO − and MOH groups. In this case, since the additional reaction is given by the combination of two elementary ones (where an exchange of ions occurs) the resulting reaction constant K c is not strictly speaking a dissociation constant and it is a dimensionless quantity.

III. SELF-CONSISTENT SIMULATION FRAMEWORK
The procedure described in the previous sections applies to a single type of site. If we assume different types of sites coexist and, from a chemical point of view, independently interact with the ionic species in the electrolyte, then they can be considered separately. The total surface charge is then the sum of the contributions by all different sites. As an example, let us consider the charge accumulated at equilibrium in M different sites. Starting from (7), the expression of the total surface charge Q S,tot can be generalized as where the probability of the site/ligand " j " to be in state "i," f i, j , generally depends on ionic concentrations and reaction constants and the M matrix stems from the block assembly of the M matrices of each site type. Note that although all sites are treated independently of a chemical point of view, their states are mutually coupled by the sensor electrostatics [15]. Thus, neglecting steric effects and assuming that the electrolyte is an infinite reservoir of ions/analytes, an iterative self-consistent solution of the Boltzmann distributions [for the mobile charged ions with concentration ρ m , (18)] and the Poisson equation [see (19) for the 1-D case] are necessary to determine the equilibrium states of the sensor system.
In (18), ψ denotes the electrostatic potential, z i the i th ion/analyte net number of elementary charges, N sp the number of ion/analytes in the bulk of the electrolyte and [ A i B ] their concentration. T is the temperature and k b is Boltzmann's constant. In (19), ε is the dielectric permittivity while ρ f is the concentration of fixed charges throughout the device. The solution of such a system can be obtained as follows. Equation (6) plugged into (17) gives the total surface charge. The latter, together with (18) (for mobile species in the electrolyte), is inserted into (19), leading to a single differential equation where the unknown is the potential profile. Its solution yields the electrostatic potential and ionic concentration profiles in the electrolyte at equilibrium and must be found numerically. In particular, at each iteration of the self-consistent Newton loop, the ions/analytes' concentration entering the terms in Fig. 8. Simulated potential profile (right) across the sketched structure (left). The surface charge arisen from a reaction at the oxide/electrolyte interface is screened by the ions in the electrolyte. We use SB model for HfO 2 oxide with parameters reported in Table I. The electrolyte consists of NaCl at 100 mM and pH = 7.1. matrix M [see (6)] need to be updated with the potential profile found in the previous iteration, to get a new estimate of the surface charge Q S,tot . The same applies for mobile charges in (18). We include steric effects by introducing a thin dielectric layer (so-called "compact layer" or "Stern Layer" from the Stern-Gouy-Chapman theory [21]) that prevents to reach unphysically large density of ions at the interfaces. In all the simulations performed in this article, the Stern Layer is modeled as a parallel plate capacitor with capacitance 20 μF/cm 2 [22]. The assumption of an infinite reservoir is valid whenever the liquid chamber is open or much thicker than the sensing layer.

IV. RESULTS
The simulation framework described in Section III [see (6) and (17)- (19)] can be easily implemented in any programming language and different sensor structures (not necessarily 1-D) can be simulated. As application examples, Fig. 8 reports a 1-D model representative of ISFET-like devices, where a solid/liquid blocking interface is studied. A null electric field is assumed at the outer boundary of the oxide layer (i.e., the MOSFET is in the flat-band condition) while the bulk of the electrolyte is contacted by an ideal reference electrode which sets the potential to zero. Fig. 8 reports the selfconsistent potential profile for a bare oxide exposed to an electrolyte where a surface charge builds up at the interface (e.g., mSB reactions with hydroxyl groups) and a rapid decay of the potential results from the high electrolyte ionic strength. The structure includes a Stern layer at the oxide/solution interface (described in Section III). In the following, we validate our approach and propose an example of model calibration from experimental data. Although we used ISFET-like devices as example, one should note that the proposed method is valid for any structure where chemical reactions with surface sites are involved. Furthermore, we focus on the effect of the charge induced by chemical reactions on the flat-band voltage of MOS systems. In this case, including the semiconductor layer is not needed. However, the proposed model can be in principle integrated with a description of the electrostatic and transport in the semiconductor region and used to analyze the effect of the induced charge on other device characteristic. We validated our model with the experimental data of a pH sensor reported in [11]. At high concentrations of chloride ions, the authors signaled a clear response of the device and developed an ad hoc model involving a single type of site (the mSB model in [11], see Section II-E). The surface charge expression is given by (16) Bulk (B) and surface (S) concentrations are linked by the Boltzmann distribution, (18). Given the physico-chemical properties of the surface (reaction constants, number of sites per unit area) the variables that matter are the bulk concentrations of the reacting ions (namely protons and chloride ions) and the background electrolyte ionic strength. Fig. 9(a) compares our model with the simulations and experiments from [11], using the same set of parameters as in [11] (see Table I). As expected from the equivalence shown in Section II-E, the results are similar, although in [11] the electrostatics was modeled with an effective double layer capacitance instead of the numerical solution of the PB equations as done here. Fig. 9(b) highlights the effect of large Cl − concentrations at high pH where experiments tend to saturate and even show a reverse trend as a function of [Cl − ]. As suggested in [23] for the NaCl salt, this may be due to ion-pair formation at the surface that neutralizes the negative charge associated with trapped chloride ions. We model this effect by adding a fifth state at the right of the diagram in Fig. 7. This state, MOHClK, that resembles the formation of KCl pairs has z 5 = 0 and is connected to the MOHCl − by a branch coefficient [K + ]/K i . In our approach, this extension of the mSB requires only the addition of one row and column to (5), whereas lengthy calculations would be required with the conventional approach. Fig. 9(b) shows that (K i = 2.6 × 10 4 M) our model can capture the behavior at pH = 12. Similar inversion of the trend at high [Cl − ] are also at lower pH values (e.g., pH = 6), but this appears to be associated with a different phenomenon that deserves additional investigation. These results confirm the ability of the general approach to provide a framework to rapidly test different physical hypotheses in the interpretation of the data.
B. Model Calibration Using Experimental Data [12] Multiple measurements at different ionic concentrations are usually necessary to calibrate our model to experimental data. While a global identification procedure through simultaneous fitting of all parameters would be possible, here we show, based on the data found in [12], that our approach allows us to follow a step-by-step procedure based on physical hypotheses. Fig. 9. Surface potential ψ 0 at the interface HfO 2 /electrolyte as a function of the Cl − activity in the bulk of the solution. The experimental data (symbols) [11] are the measured threshold voltages V th converted into surface potentials such that ψ 0 = 0 V for pH = , a Cl − = 10 −5 M. Simulations refer to (a) mSB model and (b) extended mSB with the additional state MOHClK for both our model (solid lines) and for [11] (dashed lines).
The authors reported Na + sensing with 44-mV/dec sensitivity in the physiological range of concentrations, defining the response as the differential signal of two devices, respectively, with ("active") and without ("control") a functionalization layer [see Fig. 10(a)]. The subtraction of the control response to the active one allows us to eliminate the gold intrinsic response to pH and reveal the sodium effect on the surface potential. We proceeded as follows. Starting with the bare gold (control device) response to pH and assuming that the concentration of chloride in the electrolyte was negligible, we used the SB model to find the dissociation constants K a and K b as well as the sites' density N S that gave the best match with the experimental data, as reported in Fig. 10(b). In accordance with the values proposed in [12], we assumed that the pH of zero charge (pH zc ) is equal to 7 (i.e., the product K a K b is constant) reducing the number of free parameters to two. In such conditions, K b and N S can be used to change the response slope around the pH zc and at extreme pH values. Calibration of the additional parameter K c needed by the mSB model (see Fig. 7) was possible using the experimental response of the control device to KCl salt, as the pH remains constant during the calibration while the chloride concentration is changed. Therefore, we assume that the f 1 , f 2 , and f 3 in Fig. 7 are unaffected and so the K a and, K b . The only unknown left is K c . The value of this parameter is correlated with the Cl − concentration at which the sensor shows a visible response. Fig. 10(c) shows, once again, that a good match is found using the same value as in [12]. All values of the model parameters can be found in Table I. With the same set, we found a good agreement with the experimental data for the control device response to NaCl [see Fig. 10(d)]. Both KCl and NaCl responses have been successfully modeled with the mSB model, meaning that more complex models such as the one proposed in Section IV-A were of second-order effects. The sodium-sensitive functionalization of bare gold surfaces with (neutrally charged) thiol-terminated crown ethers was shown by Wipf et al. [12] to take place only at non-oxidized gold atoms, leaving unaltered the gold response to pH. This fundamental aspect allows us to use the same mSB parameters   [12] are depicted together with their associated graph charts. For example, active sensors have two types of sites: crown ethers and oxidized states. The former selectively binds sodium cations, whereas the latter react with protons and chloride anions according to the mSB model. of the control device also when simulating the active one. This is further supported by the good agreement between the experiments and simulations of the active device to both pH and KCl [see Fig. 10(b) and (c)]. Thus, the response of the active sensor differs from the control one only for electrolytes containing Na + . In fact, the additional sites (ligands) are meant to capture sodium ions, thus compensating the surface complexation of chloride. For simplicity, we assume that the functionalization layer has a negligible thickness so that from an electrostatic point of view the surface charge is located at the gold/electrolyte interface. The sodium and ligand interaction is composed of a single chemical reaction, hence a simple Langmuir model has been used [see Fig. 1(b)] having as parameters the dissociation constant, K d , and the surface coverage of the self-assembled-monolayer, N L . The first parameter influences the concentration of sodium ions at which the ligand starts to capture a remarkable amount of them. The second parameter defines the strength of such binding in terms of surface potential change: High N L translates into less cross-sensitivity with mSB-induced surface charge. From Fig. 10(d), given the flat response observed for the active sensor with respect to NaCl salt, it is reasonable that K d roughly coincides with the corner concentration of the control device's response. At this point, N L is just tuned until the chloride ions response is balanced by the sodium one [solid line in Fig. 10(d)]. The set of parameters giving the best fit is reported in Table II.

V. DISCUSSION AND CONCLUSION
A structured methodology has been proposed to compute the surface charge due to multiple surface chemical reactions at the liquid interface of potentiometric FET sensors. The intuitive graph description of the reaction set helps to infer the mathematical expression of the surface charge at equilibrium and to quickly introduce new reactions for the purpose of testing new hypotheses in the explanation of the data. When coupled to the self-consistent solution of the PB equation, the model successfully reproduced several experimental results from the literature and shed new light on some of the observed behaviors. Although developed with reference to uniform spatial distributions of sites, the theory is also applicable to individual sites and amenable to implementation in full 3-D simulators (see [24], [25]). This would make possible to simulate non-perfectly flat sensing surfaces and/or different shapes of the electrode itself. More in general, the methodology is applicable to a broader set of potentiometric devices than the ISFETs, here used as reference structures, so far as the transduction mechanism is due to surface chemical reactions. The general framework developed in Section II-B is applicable only to static or quasi-static analyses. It is in fact based on the graph in Fig. 1(b). To include surface kinetics and extend the theory to small signal ac or noise analysis, graphs as the one in Fig. 1(a) should be used. In such cases a general matrix formulation is possible linking the d f/dt to the states' probabilities themselves and rate constants. The full development of such a theory is under investigation and, due to the lengthy mathematical derivations, will be presented in a forthcoming publication.