Modeling of an Internal Stress and Strain Distribution of an Inverted Staggered Thin-Film Transistor Based on Two-Dimensional Mass-Spring-Damper Structure

: Equipped with a two-dimensional topological structure, a group of masses, springs and dampers can be demonstrated to model the internal dynamics of a thin-film transistor (TFT). In this paper, the two-dimensional Mass-Spring-Damper (MSD) representation of an inverted staggered TFT is proposed to explore the TFT’s internal stress/strain distributions, and the stress-induced effects on TFT’s electrical characteristics. The 2D MSD model is composed of a finite but massive number of interconnected cellular units. The parameters, such as mass, stiffness, and damping ratios, of each cellular unit are approximated from constitutive equations of the composite materials, while the electrical properties of the inverted staggered TFT are characterized by utilizing an electro-mechanical coupling relation derived from the quantum mechanics. TFTs are often used in biomedical sensors/transducers attached to human skins, and, for the purpose of simulation and validation, the boundary conditions on the interface between the TFT and the human skin were modeled as a spatially distributed sinusoidal excitation with a frequency of 50 Hz, assuming the TFT thickness is more than tens of microns. The fidelity of the 2D MSD structure in the modeling of an inverted staggered TFT is verified by comparing its simulated total displacement field with thatof a finite element analysis (FEA) model. The advantages of the MSD model include a dramatic reduction in memory use by up to 60% and faster computation times that are up to 80% lower. More importantly, the MSD model is better suited than FEA to many problems in accurate tissue modeling for medical applications, for which FEA is becoming a bottleneck. This work develops a novel modeling approach, which can be extended to other types of flexible thin film transistors.

commonplace as wearables, soft robotics and medical devices become more mainstream. These developments are getting to the point that soft circuit boards and bendable electronics are soon to require higher performance modeling as the physical flexibility causes direct impacts on the electrical characteristics of the circuit-changing its behavior. As we integrate more and more complex functionality into multi-physics modelers with deeper and deeper interactions, the computational complexity and accuracy of modeling is resurfacing as a limiting bottleneck. We present a new method of modelling flexible circuit elements using MSD models as a replacement for traditional FEA. On a simple TFT, we show our model achieves equivalent accuracy, while reducing memory usage and computational time by 60% and 80%, respectively.
As a representative flexible electronic component, the thin-film transistor (TFT), has attracted much attention in its manufacturing as well as applications. TFT is mainly composed of three terminals, source, drain and gate, as well as two additional layers, namely semi-conductive and dielectric layers. In the general structure, the gate island is separated from the semi-conductive layer by dielectric insulator, while the source and drain islands are aligned on the two sides of the semiconductor. There are four commonly used TFT layouts. As shown in Fig. 1, the TFT is called staggered if the gate electrode and the source/drain electrodes are on the opposite sides of the semi-conductive layer, and the TFT is called coplanar when they are on the same side. Moreover, if the gate terminal is configured on the bottom of the dielectric insulator, the TFT is considered in an inverted layout [2]. The choice of the design is often influenced by fabrication constraints, such as the availability of chemicals and processes, as well as interaction of individual layers [3]. Among these four layouts, the inverted staggered TFT showed improved performance over other TFT layouts in the field effect mobility and the I on /I off ratio [4], and it has also been instrumental in our understanding of the multiphysics modeling and the associated manufacturing limitations for flexible electronics. Therefore, we place our interest in the characterization of mechanical and electrical properties of the inverted staggered TFTs. Multiple modeling approaches have previously been explored by researchers. For example, stochastic partial differential equation (SPDE), which is widely employed to predict the random variations of the market price for financial assets, has been applied to model the deposition process of thin films [5,6]. Several comprehensive SPDE models, such as Edwards-Wilkinson model, Kardar-Parsi-Zhang model, and ballistic deposition model, consider relaxation effects on the deposition surface, improving the precisions in predicting the printed surface roughness, the optical reflectivity, electrical conductivity and physical porosity of the thin films [7]. However, these SPDE models lose their fidelity when more important TFT's electrical properties, such as charge carrier mobility, and threshold voltage, need to be characterized. For an inverted staggered TFT, the increased thickness, the bonding effects and inhomogeneity among different layers will influence its stress distributions, and therefore will cause strained effects on its electrical properties as well. In order to precisely characterize the TFT's electrical properties, much of the current research pays particular attention to compact modeling of Silicon-based MOSFETs. For example, Sheu [8] comprehensively investigated the layout-dependent effects on Nano-scale MOSFET and proposed a stress-dependent dopant diffusion model for its numerical simulations. Wacker et al. [9] presented a DC circuit simulation using BSIM3v3 transistor model to explore variable uniaxial mechanical stress effects on CMOS transistors. A similar study in this area is the work of Alius et al. [10], who expressed the strained variations on electrical parameters as a linear function of the stress, i.e., Δ(parameter) parameter = sensitivity × stress, and used bent-chip experiments to verify this linear relationship. In a follow-up study, Heidari et al. [11] proposed a bendable MOS compact model and experimentally verified that the stress-induced sensitivity coefficients for the drain current and the threshold voltage are both equal to the modified piezoresistive coefficients of the doped Silicon. This observation contributes to many more experimentally identified statistical formulas developed for TFT's modeling research, and a detailed review work which discussed the bending induced electrical response variations on TFTs was reported by Heidari et al. [12]. Although these compact models of the Silicon-based MOSFETs accurately examine additional electrical characteristics in contrast to SPDE models, there is still very little scientific understanding of the relationship between strained variations on electrical properties and the stress distributions in the semiconductor channel of the TFT. In addition, up to now there has only been a very limited number of simulation studies of TFTs under dynamic boundary conditions. Most of previous research on TFT's electrical characterizations solely focused on the bending case. However, more complicated scenarios, such as compressing, stretching and twisting, were not discussed. It is our motivation to fill this knowledge gap by proposing a two-dimensional mass-spring-damper (MSD) structure and generalizing the boundary conditions on the MSD boundary layers to develop a holistic modeling of the inverted staggered ZnO TFT.
The MSD structure is widely used in studies for phonetics and vibration systems. However, recent research demonstrated that when the stiffness, mass, damping coefficients and the topology of the MSD structure are chosen prudently, the MSD structure can be applied for modeling soft tissues, buffered damping systems, running and hopping process of a human body [13][14][15]. It has also been demonstrated in these applications that in contrast to FEA models, MSD models exhibit better modeling accuracy, especially for coping with viscoelastic and nonlinear materials in flexible and microscale structures. Besides that, it is also reported by Chen et al. [16] that MSD model has simpler formulation, faster calculating speed and more realistic rendering performance for modeling soft tissues than FEA. This further proves the validness of MSD model in the applications for TFT modeling. In this paper, each material layer of the inverted staggered TFT is divided into a number of cellular units in repetitive patterns. The cellular unit is formed by concentrating mass to the center with four spring-damper pairs extending radially from the center with 0 • , 45 • , 90 • and 135 • respectively. The topology of the mass-spring-damper structure is determined by the quantities of the concentrated masses, the orientations of spring-damper pairs and the connection density of the mass center (i.e., the number of spring-damper pairs connecting to the same mass center divided by the geometric volume of a cellular unit). As depicted in Fig. 2, the tension and vibration on the skin surface are two sources of external excitations on the TFT's bottom layer. For a substrate with thickness around tens of microns, its adhesions to the skin surface is very small [17], which also makes the effects of surface tension negligible. In order to investigate the stress/strain distributions induced by boundary vibrations, the frequency of 50 Hz, which was provided as a typical frequency value for vibrations on human skin surface by Kandel et al. [18], is simulated in our study. Meanwhile, the electro-mechanical coupling relation, i.e., the linear relationship between strained variations on electrical parameters and the stress distributions, are fully explained from the perspective of quantum mechanics, while other typical electrical properties, such as threshold voltage (V th ) and drain-current (I D ) of the inverted staggered ZnO TFT, are characterized to reflect the simulated sinusoidally-developed excitations on the TFT and human skin interface. In a comparison study with FEM, it is revealed that the total displacement distributions in the two-dimensional MSD structure will approach the total displacements in the FEM model as the longitudinal and transversal number of cellular units increases several times over. The evaluation of MSD model highly depends on the accuracy of the simulation results from the control group, i.e., the FEA simulations. The mesh refinement study verifies the convergence of FEA simulation results as the mesh size is properly tuned. Within the scope of simulation, this observation is sufficient to demonstrate the acceptable accuracy of the FEA simulations. Accordingly, the comparison study between MSD and FEA simulations further confirms the validity of the two-dimensional MSD structure applied for modeling the mechanical deformations in the TFTs. In solid mechanics and FEA models, it is generally assumed that the underpinned materials in the structures are linearly elastic and inviscid [19]. Due to this fact, the following assumptions should be presumed: • The material is homogeneous and isotropic, and a linear constitutive relation (Hooke's law) is satisfied, the fourth-order stiffness tensor (elastic moduli) s is symmetric; • The electrical properties in the thin film are independent of temperature variations (i.e., the electrical system is decoupled from thermal system); • There is no body force (b = 0), and no heat source (ϕ = 0) in the TFT.
The two-dimensional MSD model can deal with viscoelastic, nonlinear and flexible materials by varying its damping coefficients and structural topologies, therefore we only assume that there is no body force, no heat source in the MSD model.

Mechanical Model Derived from Solid Mechanics
For an archetypical solid mechanics system (in Eulerian description), there are 23 unknown scalar variables to be determined in each spatial point: mass density ρ, displacement variables u i , i = 1, 2, 3., second-order stress tensor σ ij , i, j = 1, 2, 3., second-order strain tensor ε ij , i, j = 1, 2, 3., and temperature T. In solid mechanics, these 23 unknowns are determined by a group of interrelated partial differential equations. The tensor forms of these equations are presented in Eqs. (1)- (6). The first four equations are basic balance laws for mass, linear momentum, angular momentum and energy. Eq. (5) provides strain-displacement relationship and Eq. (6) is the materials' constitutive equation, where all these equations satisfy Einstein's indicial summation convention, c v is the constantvolume specific heat of the material, k is the material's heat conductivity, s ijkl is the fourth-order stiffness tensor, and β ij is the thermal-stress coupling coefficients in Duhamel-Neumann relation [20]. Consider the symmetry of σ , ε and s, if we let β ij = 0, Eq. (6) for homogeneous and isotropic materials can be rewritten as Eq. (7), ⎡ where E is the Young's modulus and ν is the Poisson's ratio, Eqs. (1)-(7) are the close-form of governing equations to solve all those 23 scalar variables for a continuum body. If the domain has simple geometry and special boundary conditions, the analytical solutions can be solved [21,22]. For complex geometry and boundary conditions, Eqs. (1)-(7) are solved by FEA.

Mechanical Model Derived from Mass-Spring-Damper Structure
The serial structure of 1D MSD model is shown in Fig. 3. If we take infinite 1D cellular units in this serial structure, the boundary conditions on the two ends perish. For this infinite geometrically periodic structure, Bloch theorem implies that the system states should also be geometrically periodic. The dynamics of a 1D mass-spring-damper system with infinite geometrically periodic structure can be described in Eq. (8), where u p+j is the vibration displacement for the j-th cellular unit in the p-th block, m j , c j and k j are respectively the mass, damping coefficient and spring stiffness constant for the j-th cellular unit in the p-th block. Jensen proved that Eq. (8) has an analytic sinusoidal-form solution u p+j (t) = A j e (p+j)γ −iωt , where A j is the vibration amplitude of the j-th cellular unit in the p-th block, γ is the wave number and ω is the wave frequency [23].
where f j is the external force applied to the j-th cellular unit, N is the number of cellular units in the 1D serial MSD structure. Consider the external forces have a steady vibrational frequency, i.e., f j =f j e iωt and u j =ũ j e iωt , we can write Eq. (9) in a more compact form, as given in Eq. (10), where M, C and K are N × N matrices in Eqs. (11)-(13), The mass, spring stiffness and damping coefficients for each cellular unit are predetermined by material's constitutive equation. For example, m j is determind by material's density and number of cellular units, k j is proportional to material's elastic modulus, and c j is determined by viscosity of the material. Correspondingly, the damping coefficients c j in 1D MSD model is given by Eq. (14), (14) wherek j = m j ω 2 j is the equivalent stiffness with ω 2 j = k j + k j−1 m j in 1D case, and ζ j is the damping ratio determined by material's viscosity.

Figure 4:
The structure of a two-dimensional MSD system with its inclusion filled with type-1 material and outer two layers with type-2 material (left), and the corresponding topology of a cellular unit in the two-dimensional MSD structure (right) If we analyze the two-dimensional MSD structure shown in Fig. 4, a group of second-order ODEs as given in Eqs. (15)- (16) can be obtained to describe the displacements at the discrete cellular units, where u j,k and v j,k are scalar components of the 2D displacement vector, f j,k,x and f j,k,y are scalar components of the external force applied to m j,k . In order to solve Eqs. According to the second-order elliptic PDE theory, the boundary conditions can be classified into three categories, i.e., Dirichlet condition, Neumann condition and Robin condition [21,24,25]. In this two-dimensional MSD structure, the boundary conditions can be similarly generalized to three types: D-type condition, N-type condition and R-type condition, which are respectively given in Eqs. (17)- (19), where In two-dimensional MSD structure, we let k j,k,2 = k j,k,4 = 1 2 k j,k,1 = 1 2 k j,k,3 , and the damping coefficients are determined by Eq. (20), wherek x,j,k = m j,k ω 2 x,j,k andk y,j,k = m j,k ω 2 y,j,k are equivalent stiffness in the x and y directions, and ω 2 x,j,k and ω 2 y,j,k are defined in Eq. (21). Correspondingly, we let c j,k,2 = c j,k,4 = c j,k,1 + c j,k,3 4 .

Electro-Mechanical Coupling Relation Derived from Quantum Mechanics
In order to establish electrical-mechanical coupling relation in the inverted staggered TFT, we need to borrow several notations from crystallography. First, we use the triple h, k, l , i.e., the Miller index, to denote the lattice planes in a Bravais lattice. The normal vector of plane h, k, l are the unit basis vectors of the reciprocal lattice. Since it is an industrial convention to design the CMOS devices based on the standard 001 plane, a commonly used layout for source/drain islands and semiconductor channels in the transistor is depicted in Fig. 5. In this representation, the directions of the channel currents J and the average stress σ in the plane is denoted by angle θ and ϕ, respectively.
where Π(θ, ϕ) is the piezoresistive coefficient dependent on stress/current directions in 001 plane, and it is determined by Eq. (23) [9], Π (θ, ϕ) = Π 11 · cos 2 θ · cos 2 ϕ + sin 2 θ · sin 2 ϕ + Π 12 · cos 2 θ · sin 2 ϕ + sin 2 θ · cos 2 ϕ + 2 · Π 44 · sin θ · cos θ · sin ϕ · cos ϕ where Π 11 , Π 12 and Π 44 are three fundamental piezoresistive coefficients in the semiconductor. However, as we discussed before, the neglection of stress component in the transversal direction (i.e., in (0, 0, 1) direction.) will influence the accuracy of the coupling relational model. Without loss of generality, we assume θ = ϕ (i.e., the carrier flow direction coincides with the stress component in the 001 plane.). If we further let the Manhattan style be MOSFET's manufacturing layout which requires that θ = 45 • , the piezoresistive coefficient in this case will be Π (45 • , 45 • ) = (480 ± 4) × 10 −12 Pa −1 for doped Silicon [9]. Under these settings, in our model we treat the stress σ in Eq. (22)  In quantum mechanics and solid-state physics, the threshold voltage of an n-channel CMOS device is defined in Eq. (24), where V fb is the flat-band voltage, m (∼1.2−1.4) is the body-effect coefficient, and ψ s is the potential difference between the extrinsic and intrinsic Fermi levels. The body-effect coefficient is slightly dependent on strains in the channel, nevertheless, the stains greatly influence the net shifts of the valence bands and the conduction bands for the semiconductor, leading to shifts of the extrinsic Fermi levels and the bandgaps. Quantum theory proved that the threshold voltage shift in the strained n-type MOSFET is determined by Eq. (25) [26,27], where ΔE c is the shift of the conduction band due to strains, ΔE g is the change in semiconductor energy gap, and N v is the effective density of states (DOS) in the valence bands, e is the elementary charge of an electron. The contribution from the changes in effective DOS in valence bands is very small for most semi-conductive materials, so that the last term in Eq. (25) is negligible. In addition, deformation potential theory is applied to obtain the values of ΔE c and ΔE g . The lattice scattering and external forces cause strains, which induces the potential shift. The deformation potential theory in solid-state physics shows that the band shift as shown in Eq. (26) is proportional to mechanical strains in the semiconductor [27], where Ξ are deformation potentials, and ε is the strain in the semiconductor channel. Under previously defined crystallography settings, it can be shown that the shift in the conduction band is given by Eq. (27), where1 is the unit tensor,â i is the unit wavevector (a.k.a., k vector) of valley i, { } is the dyadic product, Ξ d and Ξ u are the dilation and uniaxial deformation potentials at the conduction band edges, and ε is the strain tensor. Tab. 1 summarizes several groups of deformation potential values which have been derived and experimentally calibrated in earlier literatures. The superscripts in Tab. 1 denote (1) [27], (2) [28] and (3) [29]. Values with asterisk mean d + where a is Pikus-Bir deformation potential. All quantities in Tab. 1 are in eV. In accordance, the average shift in the conduction band is determined by Eq. (28) [27].
In order to quantify the threshold voltage shift, we also need to determine the bandgap potential ΔE g . Similar analysis reveals that ΔE g is also proportional to the strains in the semiconductor channel. Ojha et al. [30] states that ΔE g with stain is as few as meV in the advanced CMOS technology, thereby we can neglect the second term in Eq. (25). Combining Eqs. (25) and (28), the strained shift on the threshold voltage can be solved. However, it should be noticed that the deformation potentials are dependent on the semiconductor material types and the doping conditions in the channel. Moreover, the transistor layout will affect the theoretical values of the deformation potentials as well. Although the deformation potentials in the semiconductor channel are indeterministic in most real cases, there is still a significant observation according to Eqs. (25) and (28) that the stress-induced ΔV th is proportional to the average stains in the channel, so that it should be proportional to the average stress as well, as given in Eqs. (29), where λ 1 and λ 2 are two constant coefficients, which can be linearly regressed from experimental measurements. The derivations for the linear proportionality of ΔV th and mechanical stresses exhibited in Eq. (29) fully explains the linear electrical-mechanical coupling relation in previous experimental models [10][11][12]. Theoretically speaking, λ 1 and λ 2 are the sensitivity coefficients which depends on the deformation potentials in the semiconductor channel.  Considering a MOSFET in enhancement mode, the algebraic model of the MOSFET has three operational regions, i.e., subthreshold region, linear region and saturation region [31].
The drain current in subthreshold region is almost zero, and is described by Eqs. (30) and (31) in the linear region and saturation region, respectively [32], where μ eff is the effective carrier mobility, C ox is the per unit area gate oxide capacitance, W is the channel width (lateral size), L is channel length (longitudinal size), λ is the channel-length modulation parameter (∼0.01), V GS is the potential difference between gate and source and V DS is for drain and source. The strained effects on other electrical parameters are not discussed in this paper, relevant parameters (C ox , W , L, λ) are treated as constants in the simulation. Heidari et al. [11] briefly discussed strain effects on C ox in Verilog-A model, where it showed that the strained effects on C ox is infinitesimal and can be omitted. Therefore, Eqs. (22) and (28)-(31) form the computational electrical model for simulating an inverted staggered TFT under strains.

Simulation of the Inverted Staggered TFT
The structure of the two-dimensional inverted staggered TFT is shown in Fig. 6. Terminal electrodes are made of Indium-Tin-Oxide (ITO), the semiconductor is ZnO, the insulating layer is SiO2, and the encapsulating layer is polyethylene (PET). The key properties for these materials are presented in Tab. 2. The longitudinal dimension of TFT is 44 μm and the thickness is 36 μm. The thickness for each layer is also recorded in Tab. 2. The length of ZnO channel is 18 μm. For a TFT with substrate thicker than ten microns, the bonding force/adhesions on the interface between substrate and human skin surface is infinitesimal [17], and, as a result, the shear force/tensions on the substrate is also negligible. Therefore, within two sources of boundary excitations shown in Fig. 2, we only consider normal vibrational excitation in our study. In order to further strengthen our previous claim, we consider that the forearm skin can be modeled by the three-stage cascaded model in Khatyr et al. [33], where the stress-strain relation is then given by: Let us assume the channel of TFT is along the longitudinal axis of the arm. An average elastic modulus close to the longitudinal axis is reported to be E 1 = 0.657 MPa [33]. We consider E 1 as the initial Young's modulus of the skin at the elastic stage of Khatyr's model, i.e., E e ≈ E 1 = 0.657 MPa. Then it can be observed from Eq. (32) that the tension stress on the surface should be evaluated by Eq. (33), where E eff is the effective Young's modulus of the arm skin along longitudinal axis and it should be less than E e , E ve and 1/A implied by Eq. (32). Thus, we take values suggested in Khatyr et al. [33] and it is found that the tension stress on the arm skin surface is σ ≤ E e · ε ≈ 0.657 × 0.6 = 0.4 MPa. In order to observe the effects of the tensional boundary condition on the bottom layer of TFT, we compare the contour plots of longitudinal stress in the TFT with two types of boundary conditions, as shown in Fig. 7. Since the channel stress induced by the tensional boundary condition is much smaller than that induced by the normal stress boundary condition, we believe that the tensional boundary condition on the TFT-skin interface can be neglected. This is in agreement with the claim that the normal vibration is the major source of boundary excitations.
In the simulation, the inverted staggered TFT is attached on the human skin with two endpoints fixed, and a sinusoidal normal pressure (i.e., skin surface vibrations) with intensity of 20 MPa and frequency of 50 Hz is assumed on the TFT-skin interface. The choice of 20 MPa and 50 Hz in the simulation is made by referring to reported values in Kandel et al. [18] and Pawlaczyk et al. [34], with child skin elasticity modulus reported as 70 MPa and the elderly adult skin of 60 MPa. Moreover, Pawlaczyk et al. also stated that the mean ultimate skin deformation before bursting was 75% for newborns and 60% for the elderly [34]. Due to the viscoelasticity of human skin, it is implied by Eq. (33) that the intensity of normal pressure (stress) should be less than 36 MPa for elderly and 52.5 MPa for newborns. Accordingly, 20 MPa is less than 36 MPa and taken as an intermediate vibration intensity value in the simulation with acceptable and semi-empirical accuracy. The normal stress boundary condition is assumed to be a sinusoidal vibration with intensity of 20 MPa, as the normal stress on the TFT-skin interface should not be a constant value. The sinusoidal normal stress boundary condition guarantees that the intensity of interface stress can be decomposed and take values in a continuous range from 0 to 20 MPa. The intensity value 20 MPa is not the exact stress value on the skin-TFT interface. However, it should be treated as an intermediate value of the normal stress on the interface when the strain happens to be 60%. Besides that, it is found in Tab. 3 [35] that the Young's modulus of human middle back skin is 63.75 MPa in the bilinear elastic skin model (Hooke's law can be applied). If we place the TFT on the patient's middle back skin (which is a common point in electroencephalogram), 20 MPa (<63.75 × 60% = 36 MPa) is in the valid range of stress intensity for human middle back skin. Moreover, the Young's modulus of some animal's skin (e.g., cattle skin) can range from 100∼3500 MPa [36][37][38]. If we assume the TFT is placed on the animal skin, 20 MPa is also a reasonable stress intensity value for simulation study.
Even if some biometric researchers believe that the Young's modulus of human skin can be as large as 100 MPa, it is implied by Eq. (33) that the normal stress should be no greater than 60 MPa (this value is an upper bound, but not the least upper bound of the normal stress on the TFT-skin interface.). In Section 3.3, we will also discuss the stress-induced effects on TFT's electrical properties when normal stress boundary condition is applied with stress intensity of 60 MPa.

Implementation and Analysis for 2D MSD Model
The materials' constitutive equation will generate the stiffness, mass and damping coefficients for each cellular unit under different MSD sizes. In the 2D MSD model these parameters are determined as follows: m j,k = ρV /N c , where ρ is the material's density, V is the space volume of the single-material layer, and N c is the number of cellular units in the layer. Spring stiffness k j,k has four directional components. The dimensional analysis shows that their values are proportional to materials' elastic modulus. In our model we let them coincide in magnitude with elastic modulus. The damping coefficient c j,k has four components as well, their values are determined by Eqs. (20) and (21). Let M and N denote the number of cellular units in each row and in each column respectively. Then, the MSD parameters for the PET layer under different MSD sizes can be approximated and recorded in Tab. 3.
The contour plots of the total displacement in the inverted staggered TFT solved by 2D MSD model are shown in Fig. 8. In our paper, COMSOL is adopted to generate meshes for FEA, whose size is automatically adjusted by the software kernel. In order to address our viewpoint, the mesh refinement study was implemented and five groups of meshes and their FEA results are shown in Fig. 9. Additionally, the inverted staggered TFT is also simulated through FEA and its total displacement is presented in Fig. 10.
The mesh sizes and FEA simulation results processed with an 8-core Intel Core i9 CPU are summarized in Tab. 4. When mesh sizes decrease to less than 0.295 μm, the FEA simulation converges and the average stress in the active layer is solved. As shown in Tab. 4, the average longitudinal stress in the active layer is approximately 20.35 MPa. This converging tendency in FEA simulations can also be seen in Fig. 9, where the maximum field stress mainly locates within the active layer, dielectric layer and the two endpoints of the bottom boundary layer.
With the number of cellular units in the MSD model increasing by several times over, the simulated displacement field in the 2D MSD model is very close to the total displacement contours shown in FEA model. In comparison with FEA, the MSD model has much simpler formulation since the original 2D MSD dynamic equations can be compactly expressed in matrix form. With all its coefficient matrices filled sparsely, the compact MSD system reveals that less machine memory is needed in the calculations than FEA simulations. Indeed, it is found in our study that FEA simulation generally requires more than 10 GB of memory during the computations to drive the solution converge. However, it only requires less than 4 GB of memory for an MSD model with 616 × 504 cellular units to obtain a comparable solution. In addition, it should also be noted that the computational complexity of MSD is greatly reduced if the MSD systems are solved with computationally efficient numerical algorithms. Farmago et al. [39] reported that the FEA method has a computational complexity of O(NW 2 ), where W is the bandwidth of the banded stiffness matrix and N is the number of nodes, N∼1/(mesh size) 2 . This result is in agreement with the five groups of simulation times reported in Tab. 4. In contrast to FEA method, multiple iterative solvers such as GMRES, SOR and Gauss-Seidel can be employed to solve the inversion of sparse matrices in MSD model, and the simulation efficiency is tremendously improved by rendering a simulation time less than 2 minutes. Besides, the MSD approaches provide analytical expression, which can be used to derive the electromechanical relationship, as well as to perform the model reduction of the flexible circuitry, which is useful for control of electrical responses of flexible circuitry with more complicated geometrical topologies.    11 presents curves that demonstrate the linear relationship between average channel stress and the average channel displacement. The positive and negative stress indicate that the channel is under stretching and compression respectively. Under small displacements, the linear stressdisplacement relation stems from the linear strain-displacement relation in Eq. (5). Given this linear stress-displacement relation in the ZnO channel, the relationship between average channel stress and the MSD sizes can be plotted in Fig. 12. It is shown in Fig. 12 that average channel stress shows a tendency to converge as the MSD size increases. This observation also implies the validity of the 2D MSD structure for modeling the mechanical deformations in the inverted staggered TFT.

Analysis for Strained Effects on the Carrier Mobility and the Threshold Voltage
With Eqs. (22), (23) and (28), the displacement, stress and strain distributions obtained from the 2D MSD model can be utilized to predict the strained variations on the charge mobility and the threshold voltage (μ eff and V th ) in the inverted staggered TFT. Since the semiconductor in our TFT is ZnO, the piezoresistive coefficient is given as is the gauge factor of the thin-film ZnO [40]. As given in Janotti et al. [41], if the uniaxial deformed effect is neglected, the potential deformations in ZnO channel are Ξ d = −3.1 eV and Ξ u ≈0. Consider that trace (ε) = ε 1 + ε 2 + ε 3 , the strained variations on the charge mobility and the threshold voltage are displayed in Fig. 13

Electrical Characterization of the Inverted Staggered TFT
Since the TFT is assumed to be an n-type in enhancement mode, Eqs. (30) and (31) can be employed to simulate the transconductance curves of the transistor in linear and saturation regions. The strain-free parameters for our inverted staggered TFT are given in Tab. 5, where μ eff is stress-free charge mobility in cm 2 /(Vs), V th is the stress-free threshold voltage in V. The vacuum permittivity ε 0 is in m −3 s 4 A 2 /kg, ε r is the relative permittivity in ITO, t ox is the thickness of the ZnO layer in μm and C ox is the per unit area gate oxide capacitance in F/cm 2 , and C ox = ε 0 ε r t ox . Besides, we assume the gate oxide layer and the ZnO layer have the same thickness in the modeling. Accordingly, the transconductance curves of the drain current and the drain-source voltage are characterized in Fig. 14. Similarly, the transconductance curves of the drain current and the gate-source voltage are displayed in Fig. 15. As shown in these figures, the drain current I D will deviate from the strain-free value for less than 0.03%.
There are two reasons accounting for this trivial shift on the TFT's electrical characteristics. Initially, we assume that the two endpoints of the skin-TFT interface are fixed, this constraint will prevent the slipping motions between the substrate and the human skin in the longitudinal direction. Additionally, the thickness of the substrate (i.e., the PET encapsulating layer) is too large in contrast to other functional layers, so that the surface adhesion (bonding force) between skin and substrate is too small and the substrate tensions are neglected in our model. Per the analysis, it is found that a safe way to ensure that the electrical performance of the thin-film sensors is not adversely affected by mechanical vibrations or deformations is to increase the thickness of the encapsulating layer. Adding a thick encapsulating layer can not only reduce the effects of vibration, but also help move the neutral stress plane to the sensor's location if bending occurs. This finding is usually helpful for guiding the manufacturing of flexible thin film electronic devices. Table 5: Stress-free electrical parameters for ZnO channel [42] μ eff V th ε 0 ε r t ox C ox In the original analysis, 20 MPa is an intermediate value within the normal stress spectrum on the TFT-skin interface. If we assume the intensity of boundary stress takes an upper bound value (e.g., 60 MPa), the stress-induced effects on μ eff and V th will be reinforced in contrast to boundary vibration with 20 MPa. According to simulations, it is found that the relative change of μ eff is 0.1%, and 0.024% of V th . Even undergoing a boundary vibration with an upper bound intensity of 60 MPa, the enhancement-mode TFT characteristics (e.g., I D ) will be shifted only by less than 0.1%. This implies that the electrical characteristics of the TFT will not be severely impaired by the channel stress when it is attached on a patient's forearm skin.

Figure 15:
The transconductance curve between drain current and gate-source voltage for the inverted staggered TFT under a sinusoidal boundary excitation on the skin interface (t = 0.005 to t = 0.015 is a half period of the boundary excitation)

Conclusion
Prior researchers have developed multiple models to simulate TFTs. However, those models either strongly depend on parameters that need to be identified from experimental measurements, or are not capable of dealing with complex boundary conditions. This paper aims at developing a two-dimensional MSD structure to explore the mechanical deformations and electrical characteristics in an inverted staggered TFT. From the modeling and simulation of inverted staggered TFT, we theoretically proved the linear electrical-mechanical coupling relation from the quantum mechanics. Under a sinusoidal boundary excitation with intensity of 20 MPa and frequency of 50 Hz on the skin-TFT interface, it is found that the transconductance curves of the TFT shift only by 0.03% in comparison to stress-free scenarios. This finding implies that the stress-induced effects can be reduced by properly adjusting the thickness of the substrate in TFT. Meanwhile, by comparing MSD's displacement contour plots with that of a finite element model, we also confirm the validity of the 2D MSD structure in the modeling for an inverted staggered TFT and identify the advantages of MSD model over FEA model in the modeling of TFT.
These findings not only have significant implications for the understanding of how to characterize the electrical properties in a strained inverted staggered TFT, but also have a bearing on the development and innovations for modeling, design and manufacturing of general flexible electronic devices. The work scope of this paper is developing a novel modeling approach of an inverted staggered TFT attached on the human skin based on two-dimensional MSD structure. Its experimental portion for characterizing the real transconductance curves will be presented in the follow-up papers. The study conducted here certainly shed new light on the modeling research for flexible TFTs and enrich the knowledge base from the perspective of physics.
Funding Statement: This work was supported in part by the National Science Foundation through grant CNS-1726865 and by the USDA under grant 2019-67021-28990.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.