New Simulation Method for Dependency of Device Degradation on Bending Direction and Channel Length

The dependency of device degradation on bending direction and channel length is analyzed in terms of bandgap states in amorphous indium-gallium-zinc-oxide (a-IGZO) films. The strain distribution in an a-IGZO film under perpendicular and parallel bending of a device with various channel lengths is investigated by conducting a three-dimensional mechanical simulation. Based on the obtained strain distribution, new device simulation structures are suggested in which the active layer is defined as consisting of multiple regions. The different arrangements of a highly strained region and density of states is proportional to the strain account for the measurement tendency. The analysis performed using the proposed structures reveals the causes underlying the effects of different bending directions and channel lengths, which cannot be explained using the existing simulation methods in which the active layer is defined as a single region.


Introduction
Flexible displays have attracted considerable attention as next-generation displays that are flexible, light, and not easily breakable [1]. As the active layer in flexible displays, many materials have been used, such as amorphous silicon [2,3], low-temperature polycrystalline silicon [4,5], semiconductor oxides [6,7], and organic materials [8,9]. Amorphous indium-gallium-zinc-oxide (a-IGZO) films offer the advantages of high mobility, small sub-threshold swing, low leakage current, and good uniformity owing to their amorphous phase, which makes them suitable for manufacturing large-area displays [10]. Moreover, a-IGZO thin-film transistors (TFTs) exhibit superior stability against bending stress than silicon-based TFTs [11].
The degradation of a-IGZO TFTs under bending stress has been studied extensively under various bending stress conditions, namely tensile or compressive stress [6,12,13]; bending direction perpendicular or parallel bending to the current flow [14,15]; bending radius [16][17][18]; and bending cycles [19,20]. Several research groups have extracted the density of states of the active layer by conducting computer-aided design simulations and verified that the number of donor-like trap states increases with repeated bending [15,21,22].
However, the existing research has not covered the effects of device geometry on the differences in terms of the magnitude of degradation and variation of density of states (DOS). Moreover, the general method of calculating strain assumes a simple onedimensional structure and outputs a single strain value [23], irrespective of the bending direction and channel length. As can be inferred from the studies in the literature in which various bending radii have been used, a larger strain causes more severe device degradation [8,12]. Strain distribution is one of the important factors explaining this difference. Therefore, it is necessary to investigate the intensity and pattern of strain distribution in a device under each stress condition.
In this study, we conduct a mechanical simulation to obtain accurate strain distributions. Based on these distributions, a new device simulation method for flexible TFTs is suggested. The new method accounts for the dependency of device degradation on the bending direction and the channel length.

Mechanical Simulation Methods
As shown in Figure 1, a three-dimensional mechanical simulation structure was used to determine strain distribution. Multiple oxide and nitride buffer layers were placed alternately above the polyimide (PI) substrate. The bottom-gate electrode, gate insulator, and active layer were defined. An etch stopper was used to cover the a-IGZO, and the source/drain electrodes were placed on the left and right sides overlapping the active layer. Then, all these elements were passivated using an oxide film, except for the contact hole ( Figure 1b). The material properties of each layer are summarized in Table 1. tions. Based on these distributions, a new device simulation method for flexible TFT suggested. The new method accounts for the dependency of device degradation on bending direction and the channel length.

Mechanical Simulation Methods
As shown in Figure 1, a three-dimensional mechanical simulation structure was us to determine strain distribution. Multiple oxide and nitride buffer layers were placed ternately above the polyimide (PI) substrate. The bottom-gate electrode, gate insulat and active layer were defined. An etch stopper was used to cover the a-IGZO, and source/drain electrodes were placed on the left and right sides overlapping the act layer. Then, all these elements were passivated using an oxide film, except for the cont hole ( Figure 1b). The material properties of each layer are summarized in Table 1.
The channel width was set to 50 µm, and various channel lengths of 10 µm, 30 µ and 60 µm were used to derive different strain distributions depending on the chan length. The dimensions of the other elements were kept unchanged. The TFT was plac on a metal plate in the bending simulation. Two bending directions, perpendicular ( Figure 1c) and parallel (Figure 1d), were considered, where the bending axis was eit perpendicular or parallel to the current flow, respectively.  The plate was divided into several pieces in the same direction as the bending a (Figure 2a,b) with different displacements assigned to each edge. As shown in Figure when the plate is divided into six pieces, the center of the plate is fixed, and the magnitu of displacement is calculated using the distance between the center, p0, and the edges p2, and p3.  The channel width was set to 50 µm, and various channel lengths of 10 µm, 30 µm, and 60 µm were used to derive different strain distributions depending on the channel length. The dimensions of the other elements were kept unchanged. The TFT was placed on a metal plate in the bending simulation. Two bending directions, perpendicular ( Figure 1c) and parallel (Figure 1d), were considered, where the bending axis was either perpendicular or parallel to the current flow, respectively.
The plate was divided into several pieces in the same direction as the bending axis (Figure 2a,b) with different displacements assigned to each edge. As shown in Figure 2c, when the plate is divided into six pieces, the center of the plate is fixed, and the magnitude of displacement is calculated using the distance between the center, p0, and the edges p1, p2, and p3.
When the TFT is pulled along the length direction (X-axis) in the state of perpendicular bending, the displacement along the width direction (Z-axis) is set free. In addition, in parallel bending, the TFT is pulled along the width direction (Z-axis), the displacement along the length direction (X-axis) is set free. The bending radius (R) was set to 0.48 mm. In total, 580,539 nodes and 133,502 elements were used. The strain distribution on the bottom of the active layer was investigated by conducting a static simulation in the ANSYS software environment.

Strain Distribution
The normal strain distribution at the bottom of active layer under bending is shown in Figure 3. The fact that the different strain level over the layer implies that the degradation of the a-IGZO film is nonuniform. Figure 4a,b show schematics of the active layer divided into nine regions based on the strain along each bending direction. The color of an area indicates its strain intensity, and the strain intensities increase in the order of yellow, orange, and red. The detailed values are investigated along the paths cutting the plane in the length or width direction, as indicated by the red lines ( Figure 4c).  When the TFT is pulled along the length direction (X-axis) in the state of perpendicular bending, the displacement along the width direction (Z-axis) is set free. In addition, in parallel bending, the TFT is pulled along the width direction (Z-axis), the displacement along the length direction (X-axis) is set free. The bending radius (R) was set to 0.48 mm. In total, 580,539 nodes and 133,502 elements were used. The strain distribution on the bottom of the active layer was investigated by conducting a static simulation in the ANSYS software environment.

Strain Distribution
The normal strain distribution at the bottom of active layer under bending is shown in Figure 3. The fact that the different strain level over the layer implies that the degradation of the a-IGZO film is nonuniform. Figure 4a,b show schematics of the active layer divided into nine regions based on the strain along each bending direction. The color of an area indicates its strain intensity, and the strain intensities increase in the order of yellow, orange, and red. The detailed values are investigated along the paths cutting the plane in the length or width direction, as indicated by the red lines ( Figure 4c). When the TFT is pulled along the length direction (X-axis) in the state of perpendicular bending, the displacement along the width direction (Z-axis) is set free. In addition, in parallel bending, the TFT is pulled along the width direction (Z-axis), the displacement along the length direction (X-axis) is set free. The bending radius (R) was set to 0.48 mm. In total, 580,539 nodes and 133,502 elements were used. The strain distribution on the bottom of the active layer was investigated by conducting a static simulation in the ANSYS software environment.

Strain Distribution
The normal strain distribution at the bottom of active layer under bending is shown in Figure 3. The fact that the different strain level over the layer implies that the degradation of the a-IGZO film is nonuniform. Figure 4a,b show schematics of the active layer divided into nine regions based on the strain along each bending direction. The color of an area indicates its strain intensity, and the strain intensities increase in the order of yellow, orange, and red. The detailed values are investigated along the paths cutting the plane in the length or width direction, as indicated by the red lines ( Figure 4c).   Bending stress leading to a strain of more than 2.2% over 100,000 repeated cycles was considered sufficient for cracking the active layer. In one a-IGZO TFT bending experiment described in the literature, it was mentioned that cracks occurred when a strain of approximately 2.17% was applied over 4000 cycles [24]. Moreover, the direction of crack propagation differed depending on the bending direction [20], as shown in Figure 4a,b. These results suggest that both strain and cracking affect the electrical properties of a-IGZO films, and the variation pattern of DOS can differ depending on crack orientation.   The overall strain distribution pattern differs depending on the bending directi Under perpendicular bending, the strain is concentrated in the central part of the chan length (Figures 3a,c and 5a), and there is no significant difference in strain along the wid direction, as indicated by the flat curves in Figure 5b. In the case of parallel bending, strain is the highest at the center of the channel, and it decreases close to the corners (Figure 3d-f). The strain curve is parabolic along the paths that cut the plane latera ( Figure 5c) or vertically (Figure 5d).
Bending stress leading to a strain of more than 2.2% over 100,000 repeated cycles w considered sufficient for cracking the active layer. In one a-IGZO TFT bending experim described in the literature, it was mentioned that cracks occurred when a strain of appr imately 2.17% was applied over 4000 cycles [24]. Moreover, the direction of crack pro gation differed depending on the bending direction [20], as shown in Figure 4a,b. Th results suggest that both strain and cracking affect the electrical properties of a-IG films, and the variation pattern of DOS can differ depending on crack orientation.  Bending stress leading to a strain of more than 2.2% over 100,000 repeated cycles was considered sufficient for cracking the active layer. In one a-IGZO TFT bending experiment described in the literature, it was mentioned that cracks occurred when a strain of approximately 2.17% was applied over 4000 cycles [24]. Moreover, the direction of crack propagation differed depending on the bending direction [20], as shown in Figure 4a,b. These results suggest that both strain and cracking affect the electrical properties of a-IGZO films, and the variation pattern of DOS can differ depending on crack orientation.

Device Simulation Structure
As shown above Figure 3, the magnitude of damage is nonuniform over the active layer. When the device is subjected to strain, the atomic arrangement of a-IGZO changes, and trap states are generated [12]. In particular, oxygen deficiencies serve as shallow donors of free electrons in the conduction band [25]. The increased donor-like states due to the ionized oxygen vacancies induce a negative shift in the threshold voltage [22,26]. When different levels of strain are applied to the a-IGZO layer, the corresponding changes in the trap states are different. Therefore, a uniform density of states with a single region is unrealistic in the case of a strained device. Figure 6 shows the device simulation structures used to fit the transfer characteristics measured before and after 10,000 repeated bending cycles. In the case of the single-region structure (Figure 6a), the active layer is defined as a single region with a uniform density of states and is used in the device simulation before bending because of the lack of damage at that point. Two multi-region structures (Figure 6b,c), in which the active layer is defined as consisting of multiple regions with different densities of states in each region, are used for device simulation after bending.
sive region exhibits higher strain and has a greater number of donor-lik extensive regions. The transfer characteristics depending on the variatio each region are shown in Figure 7. The default curve is the simulation c to the measurements of the device of channel length 10 µm after perpen The words, 'increased' and 'decreased', in the legend means that the nu increased or decreased by 5 ×10 16 (cm −3 ) from the default concentration and donor-like states, respectively, and the other parameters are the sam default case. The variation of acceptor-like and donor-like states in the have little effect on transfer characteristic (Figure 7a,b) while the trap st sive region control the threshold voltage (Figure 7c,d). These results indic of the lower strain region is dominant in the perpendicular structure.
Second, under parallel bending, an a-IGZO film is divided into thre (Figure 6c). According to the mechanical simulation results, it should b least nine location-dependent areas along the length and width direc However, because the low strain region determines the threshold voltag flows through the high and low strain regions, as discussed in the perpen regions near the source or drain have a dominant influence on the thresh the regions in the middle. Therefore, we focused on three areas in the f the source among the nine regions and simplified the parallel multi-regi three regions, namely the extensive and intensive regions placed in para  First, the active layer subjected to perpendicular bending can be divided into the extensive, intensive, and extensive strain regions arranged in series (Figure 6b). The intensive region exhibits higher strain and has a greater number of donor-like states than the extensive regions. The transfer characteristics depending on the variation of trap states in each region are shown in Figure 7. The default curve is the simulation curve which is fit to the measurements of the device of channel length 10 µm after perpendicular bending. The words, 'increased' and 'decreased', in the legend means that the number of traps is increased or decreased by 5 × 10 16 (cm −3 ) from the default concentration for acceptor-like and donor-like states, respectively, and the other parameters are the same as those in the default case. The variation of acceptor-like and donor-like states in the intensive region have little effect on transfer characteristic (Figure 7a,b) while the trap states in the extensive region control the threshold voltage (Figure 7c,d). These results indicate that the effect of the lower strain region is dominant in the perpendicular structure.
tively. The tail state parameters and band edge intercept densities, namely NT (cm −3 /eV) and NTD 1 ×10 19 (cm −3 /eV), respectively, and the corresponding ch decay energies, namely WTA ~0.055 (eV) and WTD = 0.05 (eV), are used. The v DOS in the multi-region structure used to fit the measurements after the app bending stress is discussed in the following section. The two multi-region structures have different electrical properties owin ent arrangements of the multi-regions, as illustrated in Figure 8. The same pro multi-regions and the same density of states were used to compare the two m structures. In the perpendicular multi-region structure, the extensive region ha inant effect on the threshold voltage. By contrast, in the parallel multi-region the intensive region had the dominant effect on the threshold voltage. This threshold voltage with the structure suggests that the bending direction is impo when the same amount of strain is induced in the device. In the following s analyze the measurements using the proposed multi-region structures.  Second, under parallel bending, an a-IGZO film is divided into three regions (Figure 6c). According to the mechanical simulation results, it should be divided into at least nine location-dependent areas along the length and width direction (Figure 4b). However, because the low strain region determines the threshold voltage when a current flows through the high and low strain regions, as discussed in the perpendicular structure, regions near the source or drain have a dominant influence on the threshold voltage than the regions in the middle. Therefore, we focused on three areas in the first column near the source among the nine regions and simplified the parallel multi-region structure into three regions, namely the extensive and intensive regions placed in parallel (Figure 6c).
The two multi-region structures have different electrical properties owing to different arrangements of the multi-regions, as illustrated in Figure 8. The same proportions of multi-regions and the same density of states were used to compare the two multi-region structures. In the perpendicular multi-region structure, the extensive region had the dominant effect on the threshold voltage. By contrast, in the parallel multi-region structure, the intensive region had the dominant effect on the threshold voltage. This change in threshold voltage with the structure suggests that the bending direction is important, even when the same amount of strain is induced in the device. In the following section, we analyze the measurements using the proposed multi-region structures.
inant effect on the threshold voltage. By contrast, in the parallel multi-region struct the intensive region had the dominant effect on the threshold voltage. This chang threshold voltage with the structure suggests that the bending direction is important, e when the same amount of strain is induced in the device. In the following section analyze the measurements using the proposed multi-region structures.

Discussion
The transfer characteristics of the devices with different channel lengths before and after 10,000 bending cycles are shown in Figure 9. The threshold voltage decreased after bending, and the amount of decrease under parallel bending was higher than that under perpendicular bending. This trend can be well calibrated using the proposed multi-region structures with density of states based on the strain distribution obtained in the mechanical simulation. Because the strain level is the highest in the central region of the device with the channel length of 10 µm under perpendicular bending, the highest peak level of donor-like Gaussian states (NGD), i.e., 1 × 10 18 , is applied in the intensive region of the perpendicular multi-region structure. Moreover, the smaller NGD of 7.3 × 10 17 is applied in the intensive region of the parallel multi-region structure ( Table 2).

Discussion
The transfer characteristics of the devices with different channel lengths before and after 10,000 bending cycles are shown in Figure 9. The threshold voltage decreased after bending, and the amount of decrease under parallel bending was higher than that under perpendicular bending. This trend can be well calibrated using the proposed multi-region structures with density of states based on the strain distribution obtained in the mechanical simulation. Because the strain level is the highest in the central region of the device with the channel length of 10 µm under perpendicular bending, the highest peak level of donor-like Gaussian states (NGD), i.e., 1 ×10 18 , is applied in the intensive region of the perpendicular multi-region structure. Moreover, the smaller NGD of 7.3 ×10 17 is applied in the intensive region of the parallel multi-region structure ( Table 2).
The device simulation results indicate that the threshold voltage shift of the perpendicular multi-region structure is smaller, even though a larger NGD is used than in the case of the parallel multi-region structure. This result can be ascribed to the different arrangement of the multi-region structure. As current flows through the extensive and intensive regions arranged in series, both regions affect the current (Figure 10a), especially the extensive region. However, in the multi-region structure arranged in parallel, these regions act as independent current paths (Figure 10b), and the effect of the intensive region is not limited by the extensive region, which causes the large shift in threshold voltage.
In terms of channel length dependency, the shorter device exhibits a larger shift in threshold voltage. Under perpendicular bending, the strain near the source, which is the extensive region, increases as the channel length increases from 10 µm to 30 µm and, eventually, to 60 µm (Figure 5b). The proportional trap states are applied to the extensive region of the perpendicular multi-region structure, as summarized in Table 2, and the simulation results agree well with the measurements.    The device simulation results indicate that the threshold voltage shift of the perpendicular multi-region structure is smaller, even though a larger NGD is used than in the case of the parallel multi-region structure. This result can be ascribed to the different arrangement of the multi-region structure. As current flows through the extensive and intensive regions arranged in series, both regions affect the current (Figure 10a), especially the extensive region. However, in the multi-region structure arranged in parallel, these regions act as independent current paths (Figure 10b), and the effect of the intensive region is not limited by the extensive region, which causes the large shift in threshold voltage.   In terms of channel length dependency, the shorter device exhibits a larger shift in threshold voltage. Under perpendicular bending, the strain near the source, which is the extensive region, increases as the channel length increases from 10 µm to 30 µm and, eventually, to 60 µm (Figure 5b). The proportional trap states are applied to the extensive region of the perpendicular multi-region structure, as summarized in Table 2, and the simulation results agree well with the measurements.
Moreover, in the case of parallel bending, the density of states was applied to the intensive region based on strain (Table 3). However, the device with L = 10 µm exhibited a larger degradation despite the smaller trap states than the devices with L = 30 and 60 µm (Figure 9b). This result can be ascribed to an increase in the activated density of states with narrower Gaussian donor-like states. To investigate the effect of the width of Gaussian donor-like states (WGD), the transfer characteristics were simulated using the same peak level of donor-like states and various WGDs (Figure 11a). The donor-like states applied in the a-IGZO film (solid line in Figure 11b) were not fully activated, and the peak level of the activated states (dashed line in Figure 11b) increased because the states were distributed over a narrower energy range. Thus, the device with L = 10 µm that used smaller and narrower donor-like Gaussian states exhibited a larger shift in threshold voltage. Table 3. Density of states of the parameters of the a-IGZO layer for fitting the measurements before and after parallel bending using the single-and multi-region structures, respectively. In a comparison of the devices with L = 30 and 60 µm, the strain level at the c center differed depending on the channel length. However, the difference decreas the source (dashed line in Figure 5d). The fact that the devices with L = 30 and exhibited similar shifts in threshold voltage can be ascribed to the same levels o near the source. Thus, the same amount of NGD was applied, and it calibrated the urements well. In a comparison of the devices with L = 30 and 60 µm, the strain level at the channel center differed depending on the channel length. However, the difference decreased near the source (dashed line in Figure 5d). The fact that the devices with L = 30 and 60 µm exhibited similar shifts in threshold voltage can be ascribed to the same levels of strain near the source. Thus, the same amount of NGD was applied, and it calibrated the measurements well.

Conclusions
In this paper, we analyzed an a-IGZO TFT with various channel lengths of 10 µm, 30 µm, and 60 µm before and after it was subjected to bending stress. The strain distributions of the device under perpendicular and parallel bending were investigated by conducting a mechanical simulation. Considering the nonuniform strain distributions of the device, device simulation structures with multiple active regions were suggested for each bending direction. The different current flow mechanisms due to the arrangement of the intensively and extensively strained regions were analyzed. The proposed multi-region structures accounted for the measurement tendency whereby the device exhibited more severe degradation under parallel bending than that under perpendicular bending, even though the bending radius was the same in the two cases. Using trap states proportional to the strain, we explained the channel length dependency based on the good agreement between the simulation results and the measurements. Therefore, we expect that the proposed multi-region structures based on strain distribution can serve as a guideline in analyses of the effects of different bending conditions and device geometries.  Data Availability Statement: Data sharing is not applicable for this article.

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