General Model and Equivalent Circuit for the Chemical Noise Spectrum Associated to Surface Charge Fluctuation in Potentiometric Sensors

This paper firstly reports a general and powerful approach to evaluate the power spectral density (PSD) of the surface charge fluctuations, so-called “chemical noise”, from a generic set of reactions at the sensing surface of potentiometric sensors such as, for instance, Ion-Sensitive Field Effect Transistors (ISFETs). Starting from the master equation, the spectral noise signature of a reaction set is derived as a function of the reaction kinetic parameters and of the interface concentration of the ionic species. Secondly, we derive an equivalent surface admittance, whose thermal noise PSD produces a noise PSD equal to that of the surface charge fluctuations. We also show how to expand this surface admittance into stair-case RC networks, with a number of elementary cells equal to the number of surface reactions involved. This admittance can be included in circuit simulations coupled with a SPICE compact model of the underlying FET, to enable the physically based modelling of frequency dispersion and noise of the sensing layer when simulating the sensor and the read-out. Validation with existing models and literature results as well as new application examples are provided. The proposed methodology to compute the PSD from rate equations is amenable to use in different contexts where fluctuations are generated by random transitions between discrete states with given exchange rates.

In fact, the detectability, also known as the limit of detection (LOD), is often defined as the rms of the input referred noise (i.e. the power spectral density integrated over the measurement bandwidth) [10], [11]. Here noise comprises all electronic, (bio)chemical and diffusion/convection sources affecting the signal transduction chain. Electronic noise is well-known and often integrated in accurate models up to the circuit simulation level of abstraction. For instance, noise due to charge trapping/detrapping mechanisms in the dielectric layer [7], [12]- [14], can be included in the model of the underlying FET, the same for the dielectric polarization noise related to the imaginary part of the dielectric admittance [15], which requires suitable RC networks to describe the frequency dispersion of the gate stack impedance. Furthermore, electrolyte noise generated by thermal movement of charged species in the bulk of the electrolyte can be mimicked with thermal noise of the equivalent electrolyte resistance [6], [16], [17].
This manuscript focuses on the chemical noise generated at the sensing surface of ISFETs due to the random binding/unbinding processes. Such noise source has often been considered too low to be observable in common experimental conditions, explaining why most of the literature focuses on the noise originated by the underlying MOS device [6], [7], [11]- [15], [18], [19]. As suggested in [20] however, the use of nanoscale sensors, where ion binding has a large effect on the device characteristics, makes this contribution much more visible than in the past. Consistently, recent studies claimed successful chemical noise measurements [17], [21], [22]. In these works, the chemical noise was visible in the frequency domain spectra, shaped as Lorentzian or 1/f profiles. New models have also been developed to explain these experimental findings [16], [17], [21]- [27]. For example, in [17] the authors assumed thermal noise of the equivalent impedance representing the charged antigen protein layer bound to the sensor's surface as the relevant noise source, obtaining Lorentzianshaped spectral profiles. Other contributions were later introduced by considering mass transfer limitations [25], [28]- [30] and charge fluctuations across the Debye length [24], which result into additional and/or frequency shifted Lorentzian contributions with respect to the former model of [17]. Similar spectra have been calculated in [23], by solving numerically the mass transfer limited stochastic problem for affinitybased biosensors. Most of these approaches, however, are restricted to a single type of surface reaction. The case of multiple concurrent adsorption-desorption reactions has been considered in a broader context in the literature. For example, in [31] and [32] multiple adsorption has been considered in noise spectroscopy-based gas sensors and modeled it as the sum of Lorentzian functions weighted to fit experimental responses. Other models for multiple surface binding have been developed in [28]- [30] based on Langevin equations including mass transfer limitations of microfluidic biosensors and applied to mass fluctuations. In the context of potentiometric biosensors, Zhang et al. [21] recently proposed to interpret the chemical noise produced by the surface reactions with an equivalent electrical circuit, whose components are calculated using energy-distributed rates for the reaction kinetics.
In this paper we expand our previous work on the rms electrical noise produced by the stochastic nature of binding/unbinding events of charged analytes at the sensing surface [33], to achieve a complete spectrally-resolved theory suited to compute the noise power spectral density (PSD) for an arbitrary set of surface reactions. We show that the PSD of the surface charge fluctuation is strictly related to the kinetics of these chemical reactions, and we report a comprehensive study of the noise they generate starting from the master equation of the interface reactions. In addition, we propose a procedure to automatically generate an equivalent electrical circuit model of the interface impedance including thermal noise generators based on our theory. We then show that such model converges to the simple surface reaction models addressed so far in literature. Note that this paper focuses on the fluctuations under stationary conditions, i.e., we assume that all thermodynamical and electrochemical processes have reached local equilibrium.
The paper proceeds as follows. Section II derives the PSD for the chemical noise of the Langmuir adsorption reaction. In Section III, we generalize the approach in order to handle an arbitrary set of chemical surface reactions in a compact and simple manner. In Section IV, we employ the proposed model to derive the equivalent surface admittance due to chemical reactions. Conclusions are drawn in Section V.

II. NOISE IN LANGMUIR REACTIONS
For the sake of unambiguously setting the notation, we start from the Langmuir adsorption isotherm [34], hereafter denoted Langmuir reaction, that models the capturing of analytes, A, by specific ligands, L, tightly attached to a sensing surface. The reaction is expressed as where k f A and k b A are the forward and backward reaction constants respectively. The equilibrium dissociation constant is We assign the occupation probabilities f L and f L A as well as the net signed number of elementary charges z L and z L A , respectively, to the unbound (L) and bound (L A) states of the ligands, which are the only two possible configurations [35]. Obviously f L = 1 − f L A , so that only one occupation probability has to be computed.
In order to calculate the frequency spectrum of the surface charge density fluctuations Q S for the Langmuir reaction we start from the expression of the time-variant state probability of the bound state that depends on the adsorption and desorption rates as: At equilibrium, being d f L /dt = d f L A /dt = 0, the solution of Eq. 2 is simply given by [35] where the superscript '0' indicates equilibrium. Eq. 3 coincides with the well-known Langmuir adsorption equation [34]. In the following, we assume that the unoccupied ligands are neutral. Therefore, only the charge in state L A has to be considered. The time derivative of the state probability deviation from its equilibrium value, can be easily solved using the Laplace transform, giving where s = j ω, with j the imaginary unit and ω the angular frequency. Note that f L A is an even function of time and 2 {.} is used to include the contribution for t < 0. Fluctuations of the state probability, f L A , translate into fluctuations of the surface charge density Q S as where q is the elementary charge and W L the sensing area. The unilateral PSD of the surface charge density random fluctuations is defined as [36]: where E[.] stands for expectation, F {.} for Fourier transform on the variable τ . Since Q S is a stationary process we can set t = 0 so that the Fourier transform only applies to Q S (τ ), which is obtained combining Eqs. 5 and 6. The result, extended to N S identical sites per unit area, reads: suggests a Lorentzian spectral profile with the knee frequency equal to ν, in accordance with other models in the literature [17], [24], [25], [37]- [39].
So far the theory has been outlined for the simple case of one Langmuir reaction. In the following Section III, we extend the model to an arbitrary number of chemical reactions.

III. GENERALIZED THEORY
In this Section, we consider the general case of one or many identical sites with multiple reactions among multiple states. 1 Other extensions of this models are possible. For example, the generalization to different clusters of identical sites each responding to a different surface reaction set can be addressed as done in [35]. Conversely, the further extension to non-identical sites responding to the same surface reaction set, but with different reaction rates (situation that gives non-Lorentzian spectral profiles [21]) would require defining a discrete set of reaction constants, i.e. extending the set of sites and re-applying the model for each site. This will be however not considered in the following application examples.
Similarly to the previous case of one chemical reaction, the following study of the noise requires solving a set of rate equations (called master equation), coupled by the kinetics of the chemical reactions. As pointed out in [33], [35] if each site has N possible states, N−1 chemical reactions can be identified among them. Let us denote f i the occupation probability of a generic state and r i j the reaction connecting state i with j , as illustrated in Fig. 1. Note that, differently from our previous work on steady state analysis [35], here it is possible to have loops in the transition diagram representation.

A. Model Derivation
The out-of-equilibrium state probabilities f i are timedependent. The contribution to the master equation of the reaction r i j can be expressed as: where dots indicate null values. If more surface reactions are present, the matrix should contain four additional elements for each additional reaction. Although Eq. 9 considers a linear dependence of the matrix elements on a single ion concentration, more complex expressions can be in principle handled. Nevertheless, we restrict the following of the paper to this simple linear dependence since many examples of potentiometric sensors are based on such kind of reactions and because it will make easier the derivation of the equivalent electrical circuit in Section IV. Each reaction thus affects two lines of the matrix, with coefficients of the same magnitude and opposite sign and with forward transfer rates that are proportional to the surface concentration of the interacting ion [ A i ]. Therefore, the N×N matrix in Eq. 9, hereafter denoted as G, has rank (N−1). The f i occupation probabilities should enforce meaning that each of the N state probabilities can always be expressed as a function of the remaining ones. The system can then be reduced to (N−1)×(N−1). By denoting f the column vector of states probabilities, we define a reduced column vector f containing all the states' probabilities but one. In the following, and for the sake of simplicity, we assume that the term removed from f when creating f is always f N . Hence, considering fluctuations around equilibrium, we can rewrite Eq. 9 compactly for the new set where we use the symbol · to indicate conventional products between vectors and/or matrices whereas in the following we will denote with × the product element by element. In Eq. 11, is an (N−1)×(N−1) matrix defined as = T ·G· R, where T is the (N−1)×N transformation that eliminates the last row of G and R is the N×(N−1) transformation matrix such that f = R· f . 2 Note that since is derived from G, it will only contain analytes' concentrations as well as forward and backward reaction constants.
Analogously to Eq. 5, the solution of Eq. 11 is obtained using the Laplace transform, which is then computed along the imaginary axis j ω to yield where I is the (N−1)×(N−1) identity matrix. Once again, the factor 2 {.} is needed to account for both t < 0 and t > 0 in the transform. In order to estimate the total fluctuation of the surface charge, all the states should be taken into account. Equation 6 then becomes where z is a N-elements column vector containing the net signed number of elementary charges associated to each state, whereas f is an (N−1)-elements column vector. Therefore, as expected, Q S (t) is a scalar and thus always equal to its transpose. This observation will help writing the final expression of the PSD of the surface charge density noise in a more compact form. By inserting Eq. 12 into 7 and considering a number N S of identical sites per unit area, one gets where is a (N−1) × (N−1) matrix containing the expected values for each state of the site, The derivation of Eq. 15 is reported in Appendix I. Eqs. 14, 15 represent a closed form expression of the surface charge density noise PSD previously shown for the Langmuir reaction in Section II. When analyzing the whole sensor, the surface charge fluctuation from Eq. 14 must be considered as a forcing term inside a self-consistent solution of the coupled electrostatics, as explained in [33] and also discussed in the following of this paper.
It can be easily shown that in terms of fluctuations of vector f (i.e. f ) the vector r is cancelled out in the transformation, giving that f = R · f .
We highlight that these expressions are general and can account for an arbitrary set of surface reactions, comprising several previously published works [17], [21], [23]. The derived model thus provides the PSD of the surface charge given the surface concentration of the analytes and the forward and backward reaction constants.
B. Model Validation 1) Integrated Noise: Upon integration of Eq. 14 in the frequency domain (or angular frequency, ω), one obtains the root-mean-square (RMS) deviation that is the total noise measured when the readout circuit bandwidth is sufficiently large. Eq. 16 together with Eqs. 14, 15 matches exactly with our previously published model of the integrated noise [33], that states For simple reactions such as the Langmuir reaction described in Section II, such equality can be proven analytically, whereas for more complicated reactions we verified it numerically, obtaining perfect agreement over a broad set of chemical reactions and electrolyte compositions (not shown).
2) Comparison With Models Using Equivalent Circuits: A few authors have associated equivalent electrical circuits to surface reactions, reproducing the AC and/or noise frequency response of the sensing surface [17], [21], [40]. In order to compare results with such models, we should recall that the circuit representation of surface chemical reactions can be compactly described by an equivalent surface admittance Y S , as shown in Fig. 2, which is in parallel with the series of the Stern and diffuse layer capacitances. Consequently, the chemical noise related to surface chemical reactions is given by the respective thermal noise of the equivalent surface admittance [41]. In Fig. 2 we use the Norton equivalent circuit, with the noise current labelled i n . Thus, the noise current PSD is where k B is the Boltzmann constant and T the temperature. The surface charge, W L Q S , is obtained by integration of i n over time. Therefore, This is consistent with the derivation of Eq. 14 in Section III-A, where the surface charge fluctuation was computed at fixed potentials (meaning that in Fig. 2 no current flows in the double layer capacitance and in the gate capacitance of the FET).
Eq. 19 allows to compare the noise PSD of the surface charge density in Eq. 14 with other models based on the The chemical noise is then modelled using the thermal noise of the resistive part, with the respective Norton current noise source, i n . The contribution of the Stern and diffuse layer capacitances are represented by C st and C dl respectively. A lumped-elements description of the bulk electrolyte (C el and R el ) is also reported for completeness. equivalent circuit analysis [21], [40]. To this end, we note firstly that for the Langmuir adsorption reaction our approach is totally equivalent to that of a circuit-based model. The analytical proof of this result is reported in Appendix II.
As for other reactions, such as for instance the well known site binding (SB) of H + ions to a metal-oxide surface [35], [42], we can readily obtain the total current noise PSD, S i n i n , by considering the equivalent circuit and the thermal noise of its resistive elements. Fig. 3 demonstrates the excellent agreement between our physical model and this procedure based on the lumped element with the parameters values used in [21] (summarized in Table I). As expected, two Lorentzian spectral profiles are distinguishable, stemming from the fact that a two-reaction system was employed. Differences on the spectral profiles due to different pH values are also correctly reproduced. It is worth noting that, however, instead of implementing the whole model in [21], we focus on effective reaction constant rather than the energy resolved distributions employed in [21] to reproduce experiments. In this respect, the values reported in Table I are the effective rate constants corresponding to a measurement frequency of 10 Hz at pH = 7 (further details in [21]).
In summary, our general model for multiple reactions and multiple sites is consistent with: a) previous models for the Langmuir reaction [23], [39], [40], [43] and b) the model in [21] that treats noise due to site binding reactions as thermal noise of the equivalent circuit resistances. However, our approach summarized by Eqs. 14, 15 is much more general, as will be exemplified with a relevant case in the following section.

C. Application Example
Based on the general model described in Section III-A, we can now explore new complex but useful case studies. In fact, the generality of the model enables the noise analysis of chemical reactions with an arbitrary number of states.
In general terms, knowledge of the ionic concentrations at the very solid/electrolyte interface (see Fig. 2) is necessary Fig. 3. Comparison of surface charge density noise PSD vs. frequency calculated using Eq. 14 (black solid lines) and Zhang et al. [21] (colored dashed lines with symbols). The plots, refers to SB reactions of a TiO 2 surface at different pH values. The two models are indistinguishable. Parameters values are taken from [21] and reported in Table I. The f −2 bar is added for comparison with the Lorentzian decay, that characterizes the high-frequency response of single reactions. In this case, a tworeaction system (SB) gives the superposition of two Lorentzian spectral profiles.  [21] to calculate the surface charge noise PSD as outlined in the master equation (Eq. 9) in Section III-A. This step entails solving the electrostatics for the device structure considered (e.g. an ISFET [44]). In equilibrium conditions, the steady state electrostatic potential and the ionic concentrations are given with acceptable accuracy by the Poisson-Boltzmann (PB) differential equation [45]. Details on how to solve self-consistently the resulting non-linear numerical problem coupled with arbitrary sets of surface chemical reactions can be found in [35].
To show the potential of our model to support the analysis of complex surface chemistry, we consider self consistent solutions of an extended version of the SB model proposed by Yates and co-workers [42]. In fact, in the original model the authors beside binding/unbinding of H + ions considered  [42]. With respect to SB, the new chained reaction involves two additional states labelled MONa and MOH 2 Cl that describe the adsorption of sodium and chloride ions from the electrolyte, respectively. ion pair formations between charged SB surface groups and electrolyte ions (e.g. Na + and Cl − ). This was modeled using an additional charged plane containing the adsorbed species.
Here, we assume that adsorbed ions lie in the same plane of the surface sites. This allows us to develop an extended SB model that introduces two additional surface reactions with electrolyte ions. The result is a model that includes four reactions involving five states of hydroxyl groups (N = 5).
The complete transition diagram of the five-states reaction set is reported in Fig. 4 for sodium and chloride electrolyte ions. As can be seen, the unprotonated state M O − can either be neutralized by sodium ions or protonated by hydrogen ions. Twofold protonated states, instead, can be neutralized by loss of a proton or by chloride ions binding.
As mentioned above, we solved the PB electrostatic problem using the parameters for the HfO 2 metal oxide reported in Table II. Other simulation parameters are reported in Table III. Then, we used the calculated surface potential φ S to derive the sodium, chloride and proton concentration at the surface for specific sample compositions, that represent the steady state configuration for the surface charge noise analysis.
Since not all the reaction kinetic's parameters in the master equation are known (to the best of our knowledge) we assumed reference values for K c and K n which are kept constant for all the simulated ionic concentrations. It should be noted that, although not calibrated quantitatively, the following analysis of the sensor noise still well illustrates how the model can support complex analysis of the noise spectra and extraction of relevant noise and reaction parameters.
The results obtained with Eq. 14 for different values of k f c and k f n , and for different electrolyte concentrations (see insets) are reported in Fig. 5. a-c for a wide frequency range, well in excess of what can be experimentally achieved with a single instrument but useful to highlight all the features stemming from the complex system chemistry and to highlight where these features become visible in the noise spectra. Fig. 5 also shows the noise PSD computed for the simplest SB model that neglects reactions involving Na + and Cl − (black dotted lines). As a general rule, we observe that the noise PSD of the five-states reaction tends to redistribute the PSD towards higher frequencies with respect to the simple SB. This is a direct consequence of the modified fractional occupation of the states: less surface sites are in the M O − , M O H and M O H + 2 states and the superposition of the two Lorentzian profiles, stemming from the SB-like behaviour, gives a smaller contribution to the noise level. The high frequency excess PSD is, instead, associated to the extra-SB reactions with electrolyte ions, that generate two additional Lorentzian profiles.
Interestingly, in the examples of Fig. 5.a,b at most one additional knee with respect to the conventional SB model is observed, and its cutoff frequency appears to depend only on k f c (Fig. 5.a). The reason is that the bias point used for this analysis results in a relatively high and positive surface potential (φ S ≈ 73mV). This corresponds to an high occupancy of the positive state M O H + 2 in Fig. 4, which makes the reaction with chloride ions more favourable than the one with sodium. To confirm this interpretation we also show results for pH = 7 (Fig. 5.c)  We further note that higher values of the forward reaction constants k f generally lead to higher corner frequencies of the Lorentzian spectrum. Moreover, larger cutoff frequencies correspond to lower amplitudes of the noise PSD at such frequency points. This is in perfect agreement with our previous model [33] (see Eq. 17), stating that the integral over  Tables II and III. the whole spectrum range does not change with the reaction kinetics.
Finally, we stress that the proposed extended SB model is here used as an application example to show the capability of our model to handle complex reactions both in the electrostatic and noise analysis. Further conclusions on the simulations results require fine calibration with experimental evidence that is not part of this work.

IV. FROM CHEMICAL REACTIONS TO
ELECTRICAL CIRCUITS In this section we use our general model to derive the equivalent compact electrical circuit of the surface admittance for an arbitrary set of chemical reactions. The goal consist in expressing the surface charge fluctuations into current densities and then take the derivative with respect to the surface potential. The obtained surface admittance can be then added to the electrical components describing the sensor, as shown in Fig. 2.

A. The Total Surface Admittance
As stated in the previous sections, chemical reactions induce charge fluctuations of the binding sites, dictated by the probability of a site to be in a given complexation state f i . As proposed by Zhang et al. [21], the charge fluctuations can be expressed by means of current densities, computed as the time derivative of the surface charge density. However, this requires using a vector notation for this quantity, differently from Eq. 13. Here, each state produces its specific surface charge density, Q Si , for a total of N contributions; being N the number of states and i the state index. Each contribution can be easily related to the state probability, For the sake of simplicity, we assume that at least one state is neutral (e.g. the last one, z N = 0 and thus Q S N = 0). This assumption allows us to obtain full rank matrices in later calculations without any preliminary transformation, 3 giving: where z is the truncated version of z (without z N ) and, again, × stands for element-by-element multiplication. A N−1 elements current density vector, J S , is then obtained from Eq. 20, where the term d f (t)/dt is the same of Eq. 11 for the perturbed case and the minus sign results from the current definition (out of the interface) adopted in Fig. 2. The next step requires linearization with respect to the surface potential. From Eq. 11, this implies that the differentiation of matrix , and thus G, must be carried out. Such matrices, however, contains the analytes' concentrations whose valence charge may differ, which means having Boltzmann functions with different prefactor associated to the surface potential.
To simplify the derivation still remaining in the most general way we expand the N×N matrix G in Eq. 9 into the sum of N×N matrices containing only backward reaction constants (G 0 ), forward reaction constants, G m , and the analyte concentration A m . The two matrices G m and A m refer to the m−th species among the total M participating to the surface reactions. That is: Let us now introduce a small perturbation to the total surface charge Q S,tot , originally at its equilibrium value, due to the combined fluctuations of all the states' probabilities f and, because of the electrostatic coupling of the concentration of ions at the surface, of A m . The A m is due to fluctuations of the electrostatic potential as we will see later.
By writing these quantities as deviations from the equilibrium values J 0 S , A 0 m and f 0 and substituting them in Eq. 21, one gets where J 0 S = 0 whereas the (N−1)×N matrix T was introduced in Eq. 11. Equation 23 can be simplified bearing that G· f 0 =0 (equilibrium condition) and neglecting second order terms, i.e. assuming A m · f ≈0. This leads to In the small signal approximation, the potential variations are much smaller than the thermal voltage, φ S k B T /q, giving that the exponential dependency of the concentrations with the potential in the Boltzmann function can be linearized as: where φ S is the perturbation on the surface potential (see Fig. 2) and ζ 1 is a N×N matrix containing the charge valence of ion species A 1 (in units of elementary charges) in the same cells where matrix A 1 has non-null elements. 4 Inserting Eq. 25 into 24 and using f = R· f , one gets 4 For example, for the generic reaction r i j shown in Eq. 9, these new matrices would be · · · ζ i · · · 0 · · · . . . . . . . . .
where dotted and non-filled cells are zero and ζ i is the valence charge of the species A i .
The key point is to rewrite the first term of Eq. 26 as a function of the current density vector. This is achieved by first using the definition of surface charge density in Eq. 20. It can be shown that q N S z × · f = W L ( × ) · Q S , where Q S is the truncated vector of fluctuations of Eq. 20 and the matrix is defined as It is worth noting that the numerical implementation of the matrix may lead to undetermined/infinite elements values when more than one state is neutral. In this case, the use of symbolic calculus can avoid such situations, allowing to proceed with the calculations until subsequent simplifications eventually lead to finite numerical expressions.
By Laplace transforming Eq. 21, the relationship between surface current and surface charge density, in terms of deviation from equilibrium values becomes Hence, using the Laplace transform on Eq. 26 and Eq. 27 we obtain the expression of the admittances vector as Equation 29 contains the admittances produced by each state. Therefore, the total equivalent surface admittance Y S,tot (in units of Siemens), is given by the sum of the elements The admittance Y S,tot represents the complex fluctuations of the surface charge with the electrostatic surface potential and is represented in Fig. 2 as a bipole in parallel to the double layer capacitance. The thermal noise associated to the real part of Y S,tot accounts for the chemical noise, whereas the imaginary part represents the additional capacitive coupling between the surface and the electrolyte. High reactance values, in fact, mean larger sensitivity of the sensor.
We highlight that the proposed approach is valid for an arbitrary number of chemical reactions involving one or more than one type of site. In the following, we show that Eq. 30 is consistent with the circuit model for SB proposed in [21] and considered in Section III-B. Furthermore, in the following subsection we derive a general approach to build equivalent circuits for an arbitrary set of reactions. The respective noise can then be computed using Eqs. 18, 19.  N for c)). The ratio ΔI S,tot /Δφ S in the Laplace domain gives polynomials identical to the results of Eq. 30. These equivalent circuits are represented with a single bipole in Fig. 2 that replaces Y s . The thermal noise of the resistors produces the current noise whose Norton equivalent is depicted in Fig. 2.

B. Equivalent Circuit Extraction
The implementation of the total surface admittance of Eq. 30 in circuit simulators can be cumbersome and requires the possibility to define generic components with a specific transfer function. In this section we simplify the procedure by decomposing the total surface admittance into a set of lumped elements, namely resistances R and capacitances C. To our knowledge there has been no attempt in the literature, to develop equivalent electrical circuits of surface reactions beyond two-reactions systems (e.g. SB). Handling analytically more than two reactions is, in fact, a lengthy and tedious task without a systematic approach. Hence, we provide the derivation of a circuit model based on the general approach described in Section IV-A.
Firstly, we observe that regardless of the surface reactions involved, the order of the s polynomial at the numerator and denominator of Y S,tot in Eq. 30 is equal to N − 1.
A simple solution to build an electrical circuit that reproduces a specific spectral response is given by Cauer ladder networks [47]. Using this approach one can decompose a rational polynomial function (such Y S,tot (s)) in continuous fractions. The iterative procedure leads to the complete synthesis of an equivalent circuit by connecting elementary cells each representing one reaction. In the case of surface reactions, passive lowpass RC networks of the First Cauer Form [47] are sufficient to reproduce any driving point admittance, as illustrated in Fig. 6.a-c.
The equivalent circuit for a Langmuir reaction is given by the RC series shown in Fig. 6.a. The expressions of R and C in terms of reaction parameters are whose derivation is reported in Appendix II. Equation 31 clearly shows that the circuit time constant   Fig. 6.b whereas colored dashed lines with symbols use the equivalent circuit of Zhang et al. [21]. All the reaction parameters [21] are given in Table I. noise PSD corner frequency) only depends on the reaction kinetics and the surface concentration of analytes [25], [37], [38], highlighting the impact of their chemical parameters on the electrical equivalent. When a two-reactions system is considered (e.g. SB), a new cell is cascaded to the Langmuir-like circuit leading to the twocells network in Fig. 6.b. Further increases of the number of states will translate into additional steps in the ladder network, as shown in Fig. 6.c.
Interestingly, the equivalent circuit of the SB model [42] developed by Zhang et al. in [21] has a different topology than that of Fig. 6.b. This highlights the non-uniqueness of equivalent circuits designed to match a given spectral response. The surface admittance computed in [21] and the one generated using our method (Eq. 30 and Fig. 6.b) are compared in Fig. 7 for the datasets reported in [21] and Table I. The perfect agreement shows that the two circuits behave identically. Note that using Eqs. 18, 19 on the surface admittance derived with Eqs. 29, 30 leads exactly to the surface charge density PSD noise already shown in Fig. 3.
The major advantage of using Cauer ladders is the easiness and simplicity in extracting the component values of the circuit. Furthermore, Eq. 29 provides a structured means to determine symbolic expressions linking the R's and C's to the physical and geometrical parameters of the system. This in turn represents a powerful tool to investigate the physical origin of the noise response and possibly to steer engineering of the surface chemistry for optimum sensor performance. The ladder network of Fig. 6 can be inserted in the equivalent circuit of Fig. 2 coupled with the model of the underlying MOSFET to allow circuit level simulations of the sensor including frequency dispersion of the gate impedance and noise.

C. Equivalent Circuit for the Full Five-States SB Model
Let us now apply the method in Section IV-B to derive the equivalent circuit for the full five-states SB reaction set in Fig. 4. The focus is on verifying that the surface charge noise PSD produced by the network matches the one based on the  Tables II-III. master equation in Section III. To this end, firstly we calculated the R and C values of the four cells (see Fig. 6.c, N = 5) for different pH and electrolyte compositions. Then, we used the noise option of the circuit simulator (SPICE [48]) to calculate the PSD of the current noise (Norton equivalent as in Fig. 2, meaning that we short-circuit the networks of Fig. 6) produced by the thermal noise of the resistive components for each equivalent circuit and finally converted the result into surface charge noise with Eq. 19. Fig. 8 reports the resulting S Q S Q S spectrum for HfO 2 metal oxide with parameters as in Table II and simulation parameters as reported in Table III. In all cases the agreement between SPICE models and Eq. 14 is excellent, showing the successful extraction of the equivalent circuits and the equivalence of the two methods for the purpose of determining S Q S Q S . One should note that this corresponds to the forced surface charge density S Q S Q S , that should be applied at the device sensing surface with the electrolyte equivalent circuit and the underlying FET connected. In this way, one can ultimately determine the effect on the sensor's response allowing the integration of our theory in a circuit design flow, e.g. for the readout circuits.

V. CONCLUSION
We provided a methodology to derive the surface charge density fluctuation PSD of potentiometric sensors such as ISFETs, starting from an arbitrary set of chemical reactions. The method extends our previous model for the integrated noise [33] and reconciles, in the frame of a general theory, some existing ad-hoc models, now appearing as special cases of our more general theory. The model assumes the existence of a large number of identical non-interacting surface sites. Coupled to data on the surface reactions kinetic parameters, it can be used to better understand the noise spectral response of potentiometric sensors in a general manner.
We furthermore bridged the gap between physics-based and equivalent circuit-based chemical noise models for potentiometric sensors, showing how the latter can be consistently  THE FIRST COLUMN INDICATES THE INSTANTANEOUS VALUE OF THE  OCCUPATION, THE 2ND AND 3RD COLUMNS THE DEVIATION W.R.T.  THE EQUILIBRIUM VALUES, WHEREAS THE 4TH COLUMN IS THE  PROBABILITY TO OBTAIN THE PARTICULAR PAIR generated by turning surface charge fluctuations into noise currents.
It is worth pointing out that the model of Eqs. 14, 15 to compute the noise PSD is quite general, and applicable to problems other than surface reactions. In fact, it is a simplified version of the forward Kolmogoroff equation [49] that can be used whenever state transitions at given rates are involved (e.g. noise due to oxide traps [43], nanopores noise [27], transport through quantum dots [50], etc.).
Together with the other noise sources such as traps in the gate dielectric (included in the FET compact model) and thermal diffusion of the ions in the electrolyte (R el in Fig. 2), our model provides a complete and powerful description to interpret spectral information related to surface chemical reactions and is amenable to be incorporated in any circuit level simulation platform.

APPENDIX I CALCULATION OF THE MATRIX
The factor f 0 (1 − f 0 ) in Eq. 8 and the matrix in Eq. 15, stem from probabilistic considerations relating the average state occupation to equilibrium state probabilities.
Let us consider two arbitrary states i and j of a given site, and denote the joint occupation probability as ( f i , f j ). When i = j , at any given time the site cannot be simultaneously in both states but only in one or none of the states. Therefore, we have only three possible cases: ( f i , f j ) = (0, 0); (1, 0); (0, 1). Table IV (upper part) collect these cases, and shows the corresponding f k = f k − f 0 k with k = i, j and the joint probability.
Hence, when i = j , the expectation of f i f j is obtained by summing the product of all terms in each row of Table IV. After simple algebra, one obtains Conversely, if the two states coincide (i.e. i = j , see bottom part of Table IV) the site is or is not in a particular state, giving only two possible configurations: ( f i , f j ) = ( f i , f i ) = (0, 0); (1, 1). Using the same procedure described above, the expectation becomes The matrix in Eq. 15 is just the matrix representation of all these cases. Therefore, expressions as in Eq. 33 are expected in the diagonal elements of , whereas expressions as in Eq. 32 are those of the off-diagonal elements. It can be easily shown that Eq. 15 is just a compact form for expressing .
APPENDIX II NOISE IN LANGMUIR REACTION AND DERIVATION OF THE EQUIVALENT CIRCUIT By using Eq. 24 for the Langmuir reaction introduced in Section II, we obtain: If only one state contributes to the surface charge (i.e. the ligand L is neutral) we can ascribe to this state the total net surface charge density, Q S (t) = q N S z L A f L A (t).
On the other hand, using the Laplace transform we have that Q S (s)W L = −(W L/s) J S (s) = −(1/s) I S (s). Hence, substituting the first term of Eq. 34 and rearranging, one gets the total admittance The equivalent circuit of the admittance in Eq. 35 is shown in Fig. 6.a. In fact, the total noiseless admittance of such a circuit is Y RC (s) = 1 R + 1/(sC) = 1/R 1 + 1/(s RC) .
By comparing Eq. 35 and 36, we obtain Eq. 31 in the main text.
For this elementary case, the equivalent capacitance, C, can also be found using the definition C = W L ∂ Q S /∂φ S , substituting the state probability f with its equilibrium value and using the small signal approximation ∂[ A]/∂φ S ≈ −qζ A [ A]/k B T , where ζ A = z L A since z L = 0.
A final consideration concerns the noise produced at the interface. As stated in Section III-B, the chemical noise is equivalent to the thermal noise of the circuit representation (see Eq. 18) due to the resistive components, namely, R.
The noise current of the admittance, i n , is then given by the current divider, where we substituted s = j ω and where S i R i R = 4k B T /R. Since the total surface charge density is obtained after integration of the surface current, Q S (s) = 1/(sW L) i n (s), we can derive from Eqs. 19 and 37 the surface charge density PSD, where we used k f A [ A] = f 0 L A ν. Equation 38 is identical to Eq. 8 demonstrating that considering the equivalent circuit and its associated thermal noise are just alternative ways to compute the chemical noise. A similar result was also derived for charge fluctuations at the channeldielectric interface in MOS structures [51].