Topology optimization for piezoresistive nanomechanical surface stress sensors in anisotropic 〈111〉 orientations

Microelectromechanical systems (MEMS)-based piezoresistive nanomechanical sensors are compact sensing platforms widely employed in vapor sensing, environmental monitoring, and biosensing. Despite their extensive utility, their lower sensitivity relative to their optical readout counterparts has been a limiting factor, constraining the wider application of this technology. Prior research has suggested that alternative silicon orientations, such as 〈111〉 orientations in (110) wafers, can significantly improve the sensitivity of piezoresistive sensors. However, the complexity of optimizing two-dimensional stress distribution and handling anisotropic elasticity has made device design a formidable task, leaving this promising avenue largely unexplored. To address this challenge, we employ density-based topology optimization to generate a series of optimized designs for piezoresistive nanomechanical sensors manufactured along 〈111〉 orientations. The properties of the immobilization layer—the functional coating on the sensor—are parametrically varied to explore optimal designs. Our study reveals a transition in optimized designs from a double-cantilever configuration to a suspended platform configuration, dictated by the stiffness ratio between the immobilization layer and the silicon layer. This transition is attributed to the shift in the neutral plane and the prevailing stress relaxation mechanism. In addition, we scrutinize the effects of piezoresistor geometry and find that the optimized designs depend asymmetrically on the piezoresistor position, a characteristic stemming from the anisotropic elasticity in 〈111〉 orientations. These optimized designs, verified by finite element analysis (FEA), demonstrate a notable improvement in sensitivity of more than 20% when benchmarked against traditional rectangular designs and equivalent optimized designs in conventional orientations, thereby validating the effectiveness of the present model. This study provides crucial knowledge for the design of piezoresistive biosensors, facilitating more efficient geometric design in future sensor development.


Introduction
The advancement of nanotechnology has catalyzed the development of various microelectromechanical systems (MEMS)-based sensors adept at detecting chemicals and biomolecules with exceptional sensitivity and specificity. Among different types of MEMS-based sensors, nanomechanical sensors have garnered significant attention due to their availability and high sensitivity [1][2][3]. Compared to the conventional optical readout, piezoresistive readout-based nanomechanical sensors demonstrate promising capabilities in numerous biosensing applications, as they eliminate the need for laser alignment and function in opaque liquid environments, facilitating in situ detection of DNA [4], RNA [5], and proteins [6,7]. Nevertheless, despite their compactness, piezoresistive readout suffers from a low signal-to-noise ratio attributable to limited sensitivity and elevated electronic noise. Prior studies have proposed various strategies to mitigate electronic noise, including miniaturization [8] and optimizing boron doping [9,10], yet low sensitivity continues to obstruct the broader application of this technology.
The sensitivity of piezoresistive nanomechanical sensors can be improved through material design. For example, substituting silicon with softer materials such as SU-8 can reduce stiffness and consequently enhance sensitivity [11,12]. However, these soft materials are prone to age effects, and their performance is heavily influenced by environmental conditions [13], rendering them less robust compared to conventional materials, such as silicon.
In the context of silicon-based piezoresistive sensors, one way to enhance the sensitivity is by using alternative silicon orientations with higher piezoresistance coefficients. It has been recommended in previous literature that silicon in 〈111〉 orientations is an effective alternative to conventional 〈110〉 orientations for piezoresistive diaphragm pressure sensors [14]. However, the elevated piezoresistance coefficients concurrently introduce anisotropic elasticity, complicating sensor design and hindering the identification of the optimal geometry that best capitalizes on high sensing performance.
Designing a sensor geometry that accommodates anisotropic material properties calls for a strategy that surpasses traditional optimization based on dimensional reduction, such as beam theory [15][16][17][18] and plate theory [19][20][21]. The response of a piezoresistive microcantilever to surface stress loading is dictated by the tensor product of stress and piezoresistance averaged at the piezoresistor. Given that the piezoresistance coefficients of a material can be either positive or negative, maximizing sensitivity entails stress maximization in some directions and minimization in others-a task too complex for simple models with a limited number of degrees of freedom. Consequently, numerical studies utilizing finite element analysis (FEA) have been conducted to establish optimal principles for cantilevers with elementary geometries, such as rectangles and triangles [22][23][24][25]. However, these studies have focused on cantilevers with rudimentary shapes and isotropic materials, leaving a vast design space involving complex boundaries and asymmetric designs uncharted. A systematic design approach is therefore essential to identify sensor geometries that exhibit the highest sensitivity and optimally exploit anisotropic orientations.
Density-based topology optimization is a computer-aided inverse design method employed across various engineering domains to address intricate optimization challenges. By discretizing a design into thousands of variables and resolving the multi-dimensional optimization problem, this algorithm can yield non-intuitive optimal designs that maximize or minimize specific objectives [26]. Pederson [27] employed this approach to generate a series of optimized designs under diverse conditions, such as varying sizes and placements of the piezoresistor and anisotropic silicon orientations. However, the optimized designs produced at the time suffered from a numerical instability known as the one-node connection problem, wherein materials in the design were linked by elements with vanishingly low density. This numerical defect severely compromised the manufacturability of the optimized design, rendering it challenging to verify optimization results using FEA and experiments, and casting doubt on the effectiveness of the optimization algorithm.
In this study, we investigate the optimal sensor design for piezoresistive biosensors fabricated along anisotropic 〈111〉 orientations. By integrating topology optimization with robust formulation, a numerical technique developed recently to efficiently address the one-node connection problem by imposing a robust length scale on the optimization problem, we successfully generate a series of optimized designs under various material properties and geometric constraints. Specifically, we modify the material properties of the sensor to explore the optimal design in different application scenarios, such as vapor sensing and biosensing, and identify the most efficient designs under anisotropic silicon orientations. Further, we vary the shape and position of the piezoresistor to examine their combined impacts on the optimized designs. All optimized designs are verified by FEA, and they exhibit markedly improved sensitivity and distinct stress behaviors when benchmarked against conventional rectangular designs under identical conditions. To comprehend the significance of various geometrical features in the optimized design, we develop an FEA model to scrutinize the origin of high sensitivity and the role of anisotropy, offering insights into the rational design of piezoresistive nanomechanical sensors. Through a comprehensive investigation of designs based on anisotropic materials, this study provides crucial knowledge for the design of piezoresistive biosensors, enabling more efficient geometric design in future sensor development.

Piezoresistive sensor
The bending of the piezoresistive microcantilever is modeled in 3D as a bimorph consisting of an immobilization layer and a silicon structural layer under a fixed-free boundary condition. The cross-section along the x-z plane is shown in figure 1. The surface stress in the immobilization layer is caused by molecular interactions, and it is modeled by a biaxial stress b s (in units of N·m -2 ) relating to the surface stress (in units of with t i being the thickness of the immobilization layer. The bending deformation caused by the surface stress is calculated by solving: where F s denotes the force vector containing the bending moment load induced by the surface stress, K is the design-dependent stiffness matrix, and u s is the deformation vector under surface stress. The constitutive relation yields the stress in the cantilever: where C is the stiffness matrix of the material [28]. The piezoresistor is modeled at the clamped end and the surface of the silicon, and it is assumed to have zero thickness to simplify the modeling. The resistance change of the piezoresistor with area A is determined by the averaged tensor product between stress and piezoresistance under unit surface stress loading, which can be simplified into the following form by assuming plane stress conditions [29]: p l s and t s are the piezoresistance coefficients and stresses along the longitudinal and transverse directions, respectively. The values of piezoresistance coefficients in these two directions depend on the crystallographic axes of the piezoresistor in the silicon substrate, and they are computed by rotating the piezoresistance tensor with respect to the standard 〈100〉 directions: where l, m, and n are known as direction cosines and they are given by the Euler angles ( , f , q y) of the corresponding wafer plane and orientation. , 11 p , 12 p and 44 p are piezoresistance coefficients of p-type silicon in 〈100〉 directions according to the Voigt matrix notation as listed in table 1. Since we model the piezoresistive microcantilever in a way that the longitudinal direction aligns with the x-axis and the transverse direction aligns with the y-axis, l s corresponds to xx s while t s corresponds to . yy s Therefore, the expression for the sensitivity of the piezoresistor (i.e., the resistance change) becomes When the piezoresistor is orientated along 〈110〉 directions that are conventionally used for piezoresistive sensing, neglecting small piezoresistance coefficients (i.e., 11 p and 12 p ), equations (4) and (5) lead to a simplified form of the sensitivity in equation (6) as where the sensitivity is maximized when the difference between stresses in the x and y directions is maximized [30]. In these directions, explicit calculation of piezoresistance coefficients using Euler angles of (0,0, /4 p ) yields 71. 8   . Therefore, the sensor output can no longer be simply treated as the difference between stresses in the two directions, but a weighted difference defined by the piezoresistance coefficients.
The previously used Euler angles are used to construct a 6 by 6 rotation matrix a to rotate the stiffness matrix to account for the anisotropic elasticity of silicon at 〈111〉 orientations according to [31]: where the superscript c denotes the stiffness matrix elements in the crystallographic coordinate as listed in table 1.

Topology optimization
Topology optimization is typically formulated in a predetermined design domain Ω with a given set of loading and boundary conditions, in which the algorithm finds the material distribution that optimizes a given objective function f [34,35]. Specifically, the design domain is discretized into finite elements where the material distribution is described by an element-wise constant density vector . j r The density of each element can take values between 1 and 0, representing the presence and absence of material, respectively. Intermediate densities are allowed such that the optimization variable is continuous. The material properties and the biaxial stress representing surface stress are weighted by the element density raised to the power of p, which is a penalization factor based on the method of Simple Isotropic Material with Penalization (SIMP) [36,37] that guarantees the algorithm to generate discrete designs with only 0 and 1 densities. p is typically 3 depending on the Poisson's ratio of the material [38].
The optimization is carried out through the following steps: (1) Evaluate the objective function of the design using FEA; (2) Evaluate the sensitivity of objectives and constraints with respect to the density variation of each element based on FEA and associated adjoint solutions; (3) Update the densities of all elements based on the sensitivity analysis using the method of moving asymptotes (MMA); (4) Repeat steps (1)-(3) until a convergence criterion is met or a maximum number of iterations is reached [26,39].

Filters and projections
Density filtering based on a Helmholtz filter [40] is used to prevent the well-known numerical instability of checkerboarding [41] and mesh dependence [42]. It is implemented by calculating the elemental density as a weighted average of all the neighboring elements within a radius of r min at each iteration: where j r and j  r are the original and filtered density of element j, respectively.
Since the density filter introduces transition zones between solid and void regions that are not physically meaningful, a Heaviside projection function is further implemented to eliminate intermediate densities [43,44]. The projected density is given by [45]  where ̅ j  r represents the projected density at element j, b represents the projection slope, and h represents the projection threshold, which determines the density at which the projection function applies.

Robust formulation
To enforce a minimum feature size for solid and void regions, robust formulation evaluates the objective functions of dilated, intermediate, and eroded designs during each iteration, which correspond to The MMA solver will automatically optimize the performance of the worst-performing design among the three. By solving this minimax problem, a minimum length scale is enforced to the optimized design, which effectively eliminates one-node connection problems and improves the robustness of the design to fabrication errors [46][47][48].

Optimization problem
The design domain W is defined as the square interface between the silicon and immobilization layer as illustrated in figure 2(a), where the interface, top, and bottom surfaces of the bending bimorph are delineated by solid and dashed lines. During the optimization process, the material distribution in the design domain evolves to maximize the objective function as illustrated in figure 2(b). The surface stress and material properties of the design are weighted by the projected density ̅ j  r via a general extrusion function, which means that the presence of an element simultaneously represents the presence of the immobilization layer and silicon with predefined thicknesses as well as local surface stress. With this treatment, the thicknesses of the immobilization layer and the silicon are assumed to be constant everywhere in the design domain. The piezoresistor is regarded as a passive area with a density of 1. The optimization objective is the sensitivity of the sensor, that is, the fractional change in resistance at the piezoresistor divided by unit surface stress : The optimization problem can be described by the following minimax problem:

Computational implementation
The topology optimization is implemented in the commercial FEA software COMSOL Multiphysics 6.1. The structural deformation of the sensor is simulated using the structural mechanics module, wherein a coupled interface of shell and solid is employed to model the stressed immobilization layer and the silicon structural layer, with the immobilization layer modeled by the shell and the silicon layer modeled by the solid. The model is discretized by linear quadrilateral elements with 50 elements per side in the in-plane directions and three elements in the thickness direction to prevent shear locking. The optimization is performed with the optimization module using the globally convergent version of MMA. The adjoint problem is solved automatically in COMSOL to obtain the sensitivities for the objective and constraint functions. COMSOL is interfaced with MATLAB to automate the parameter sweep. To ensure the optimization reaches the global optimum, a continuation method involving ramping up the projection slope b from 1 to 16 by doubling every 20 iterations is employed [39,46,50]. For the robust formulation, the projection threshold h is set to 0.3, 0.5, and 0.7 for dilated, intermediate, and eroded designs, respectively. The initial volume and the volume constraint are set to 0.6. The intermediate design is considered the final design. Other model parameters remain constant during the optimization, as shown in table 2.
To evaluate the sensing performance and stress distribution in the optimized designs, the optimized solutions from topology optimization are reimported, meshed, and evaluated in COMSOL using two times denser quadratic triangular elements compared to the optimization model to ensure accuracy. The results obtained from these verification simulations are used to analyze the sensor performance and stress distribution.

Results and discussion
To examine the optimal sensor designs made along 〈111〉 silicon under various application scenarios and design constraints, the topology optimization model is employed to generate optimized designs under various immobilization layer materials, shapes, and positions of the piezoresistor. The default values of these parameters are listed in table 3 unless explicitly specified.
To illustrate a typical optimization result, an optimized design generated by default settings is shown in figure 3, with the optimization history compiled in figure S1 of the supplementary materials. The small square at the left boundary represents the piezoresistor. The optimized design features an asymmetric 'double-cantilever' configuration, where the piezoresistor is located at the junction between the two combined cantilevers. The asymmetry is in contrast to the symmetric designs obtained in previous reports where isotropic orientations were considered [51]. The double-cantilever design leads to a concentrated stress profile at the piezoresistor, as shown in figure 3(b), the averaged value of which, according to equation (6), constitutes the signal output of the sensor. The concentrated stress profile stems from a synergy of high stress in the x-direction and low stress in the y-direction at the piezoresistor, as demonstrated in figures 3(c)-(d). The synergized stress distribution, and consequently the high sensitivity, is a result of the non-intuitive arrangement of materials enabled by topology optimization.

Material of the immobilization layer
As an essential component of a biosensor, the immobilization layer provides the chemical specificity and the driving force for the transduction of chemical signals. Different applications require distinct types of immobilization materials, and the geometry and elastic properties of the material govern the design of the sensor. To better understand the impact of immobilization materials on the design of piezoresistive sensors and enable application-specific sensor designs, we systematically varied Young's modulus and thickness of the immobilization layer materials and compared their optimized designs obtained through topology optimization. Young's modulus E i and thickness t i of the immobilization layer are varied from 10 to 1000 GPa and from 0.01 to 3 μm, respectively. This extensive range of parameters encompasses practically relevant immobilization layer materials, such as gold (with Young's modulus of 80 GPa and a thickness in the tens of nanometers) and polymeric materials (with Young's modulus of 10 GPa and a thickness of several micrometers) [4][5][6][7]52]. We omitted lower Young's moduli from our study, as optimization results exhibited minimal alteration at such low values.
The optimization results underscore a potent dependence of optimized piezoresistive sensor designs on the geometry and elastic properties of the immobilization layer, as exemplified by a selected set of optimized designs in figure 4. The full gamut of optimized designs is demonstrated in figure S2-S5. When the immobilization layer possesses a relatively low Young's modulus and thickness, the double-cantilever configuration appears to be optimal, whereas a suspended platform configuration is optimal for higher Young's modulus and thickness, as displayed in the upper left and lower right corners of figure 4. The origin of this transition in optimal configurations is not fully understood. It is clear, however, that this transition is associated with the increased rigidity of the immobilization layer and the shifted neutral plane of the bilayer structure of the immobilization layer and silicon, as implied by the diminishing bending deformation and piezoresistive signal of the sensor (figure S6).
To investigate this phenomenon, we plotted averaged bending displacement of the piezoresistor and two other key attributes of optimized designs against the stiffness ratio (r) between the immobilization layer and silicon as shown in figure 5. The stiffness ratio is defined as where E s and t s are the Young's modulus and thickness of the silicon, respectively. The product of E and t represents the in-plane stiffness of a layer. This concept comes from the calculation of uniform elongation/contraction strain c of a bimorph, as opposed to the thickness-dependent bending strain, under surface stress loading [53]: where  represents the biaxial strain of the immobilization layer. Equation (13) relates the uniform strain of the entire bimorph with the strain in the immobilization layer, with a proportional constant being the ratio of the inplane stiffness of the immobilization layer and the total stiffness of the bimorph. By substituting in the definition of stiffness ratio r, the equation (13) can be written as: where ( ) / r r 1 + is the transformation ratio that scales from 0 to 1 asymptotically as r increases from 0 to infinity. When this ratio is small, only a small portion of the strain transforms into the uniform strain. When the ratio is 1, all strain transforms into uniform strain.   The neutral plane is the plane that experiences zero strain during bending deformation, and its position t n relative to the bilayer interface is calculated using [53].
We take Young's modulus of silicon E s to be 175.1 GPa, which is the average elastic moduli of the first three diagonal elements in the stiffness matrix of silicon in 〈111〉 orientations. The results presented in figure 5 demonstrate that the optimal design of piezoresistive sensors strongly depends on the stiffness ratio between the immobilization layer and silicon. Specifically, the double-cantilever configuration is optimal when the immobilization layer has a negligible stiffness, that is, r = 1. At this point, the neutral plane is located below the immobilization layer/silicon interface and at the mid-plane of the silicon layer, with a position t n of approximately -1.5 μm. Moreover, the transformation ratio ( ) / r r 1 + is close to 0, meaning there is barely any uniform deformation. Since there are only two possible mechanisms, uniform elongation/contraction and bending [16], that can relax the strain mismatch between two layers and reach equilibrium, a near-zero transformation ratio then means that the stress relaxation is dominated by bending deformation. As the stiffness ratio r increases and approaches 0.1, the immobilization layer can no longer be ignored. This leads to a rapidly diminishing bending displacement, a rising ( ) / r r 1 , + an up-shifting neutral plane, and the transition in optimal configuration: the gap between the two cantilevers deforms and develops an L-shape configuration.
As the stiffness ratio r continues to increase and exceeds 1, the offset between the bilayer interface and neutral plane vanishes and then becomes positive, indicating that the neutral plane now resides within the immobilization layer. Moreover, the transformation ratio is now getting closer to 1, meaning most of the strain in the immobilization layer becomes uniform strain in the bimorph. The dominating stress relaxation mechanism shifts from bending deformation to uniform elongation/contraction, which accompanies the emergence of the suspended platform configuration. Therefore, the drastic transition of optimal configuration could be a result of the shifting stress relaxation mechanism as the suspended platform configuration may be more efficient in generating piezoresistive signals than the double-cantilever in the regime where stretching deformation, instead of bending deformation, dominates.
Since practically important immobilization layers made from gold and polymeric materials have a typical bending stiffness E t i i of the order of ∼10 3 and ∼10 4 , respectively, while the bending stiffness of the silicon E t s s considered here is ∼10 5 . This leads to a stiffness ratio below 0.1, making the double-cantilever configuration optimal in these cases. Therefore, we will omit the L-shape gap and suspended platform configurations in the following discussion. However, it is important to note that when the thickness of the silicon is reduced to nanometer length scales where the stiffness of the immobilization layer becomes significant, further consideration on these designs may be required as they could potentially provide better sensing performance than the double-cantilever design.
In the present study, geometric nonlinearity has been omitted due to its associated high computational cost. The optimization process involves only linear analysis; as a result, variations in the magnitudes of surface stress ( s s ) and biaxial stress ( b s ) do not influence the resulting optimized designs. Consequently, any adjustments to b s to maintain a constant s s when varying t i will not affect the optimized designs, indicating that transitions of optimized designs are purely geometric. It is important to note that nonlinearity often arises from the clamped boundary condition under compressive surface stress loading, including the fixed-free boundary condition considered in this study [54,55]. However, this omission of nonlinearity here can be justified by findings by Buhl et al [56], who reported only marginal improvements in designs upon incorporating geometric nonlinearities compared to those generated using linear analysis. Nonetheless, in scenarios where the applied surface stress is large, the silicon layer is thin, or the problem is formulated under more stringent boundary conditions, such as a fully-fixed boundary condition [57], nonlinearity will become significant enough to alter the structure's stiffness and even lead to elastic instability [58]. In such instances, nonlinear analysis becomes essential, as it is known that snap-through motions during instability can result in substantial changes in the optimization outcomes [59].

The shape of piezoresistor
In order to design a piezoresistive sensor with a high signal-to-noise ratio, it is essential to optimize the doping profile, placement, and geometry of the piezoresistor, in addition to a sensitive immobilization layer as discussed in the previous section. In this section, we scrutinize the impact of the piezoresistor shape on the optimized design and explore non-intuitive design schemes that may achieve better performances. We first study the effect of the piezoresistor shape on the optimized design by varying its aspect ratio between 1/9 and 9, with the optimization results depicted in figure 6. The stress distribution in the x and y directions is illustrated in figures S7-S9. The total area is kept constant to maintain a consistent signal-to-noise ratio, as the noise level is contingent upon the doping area [9,10]. For low aspect ratios (from 1/9 to 1/7), in contrast to the double-cantilever demonstrated in figure 3, a triple-cantilever configuration emerges. As the aspect ratio increases, the central cantilever disappears, and the optimized design reverts to a double-cantilever configuration. The disappearance of the central cantilever can be rationalized by the fact that while the double-cantilever configuration is an efficient mechanism for generating high-sensitivity regions, it has a limited area of influence, necessitating the low aspect ratio design to produce two such configurations to compensate for its widened sensing area.
Interestingly, when the aspect ratio exceeds 4, an oblique hollow feature begins to develop from the clamped end. This feature does not enhance sensitivity; instead, it emerges as an artifact caused by the volume constraint. As the aspect ratio increases, the cantilevers need to elongate to maintain the concave double-cantilever geometry. However, this growth is hindered by the scarcity of materials, i.e., the imposed volume constraint. As a result, these designs must hollow out themselves to conserve materials for cantilever extension, leading to degradation in design quality at high piezoresistor aspect ratios.
To validate the efficacy of topology optimization, the sensing performances of optimized designs are benchmarked against conventional rectangular cantilevers via FEA modeling under the same piezoresistor shape. The conventional rectangular cantilevers are modeled with a similar fixed-free boundary condition, and their length and width are adjusted to maximize the piezoresistive sensitivity according to the principles suggested by previous studies [22][23][24][25]. Specifically, the cantilever length should be minimized and set to one mesh longer than the piezoresistor to avoid the high-stress region near the free edge. Secondly, the cantilever width should be maximized and set to be the same as the design domain, as depicted in figure 7. The width of the design domain and the mesh density mirror those of the topology optimization to facilitate a fair comparison between the two design approaches. The simulation results, including the distribution of piezoresistive signals and stresses in the x and y directions, are displayed in figures S10-S12.
The sensitivities and associated stresses of optimized designs and rectangular benchmarks are assessed and presented in figure 8, which illustrates the pronounced dependency of sensitivity on the piezoresistor shape  within both design strategies. As the aspect ratio increases, the sensitivity of the optimized designs initially ascends and subsequently descends, attaining a maximum sensitivity of 13.7 10 4 -(N·m -1 ) -1 at an aspect ratio of 1/2.5 ( figure 8(a)). This is 22.3% higher than the maximum sensitivity of rectangular benchmarks, which is 11.2 10 4 -(N·m -1 ) -1 at an aspect ratio of 1/7 ( figure 8(b)).
Elevated sensitivity should result from the maximization of xx s and the minimization of , yy s as underscored in section 2.1. Indeed, in optimized designs, the average stress in the x-direction xx ave , s peaks at an aspect ratio of 1/3, while yy ave , s reaches its minimum at an aspect ratio of 2. The interplay between these two stresses yields the optimal aspect ratio of 1/2.5. Conversely, the sensitivities of rectangular benchmarks exhibit a plateaued sensitivity at low aspect ratios and consistently trail those of the optimized designs. The sensitivity of the rectangular benchmark is a monotonically decreasing function of the piezoresistor aspect ratio, signifying a deficiency of design flexibility in this conventional model.

Piezoresistor position
As illustrated in the preceding section, the optimal aspect ratio of the piezoresistor for optimized designs lies around 1/2.5, driven by the extremization of stress along the x and y directions. However, the asymmetry evident in the designs obtained thus far necessitates further scrutiny, as they suggest that materials in the upper and lower halves of the design domain may contribute to sensitivity with varying efficacies. Therefore, displacing the piezoresistor from the center might modify the optimal aspect ratio and engender alternative optimal designs.
To explore alternative designs with off-center piezoresistors, we adjust the position of the piezoresistor from the lower to the upper of the design domain. We define the offset of a piezoresistor as the distance between the midpoints of the piezoresistor and the clamped boundary. Thus, with a design domain having an edge length of 300 μm, a square piezoresistor with an edge length of 40 μm situated at the lower edge equates to an offset of −130 μm, while the upper edge corresponds to 130 μm. We investigate various positions with offsets ranging from −110 to 110 μm while the piezoresistor aspect ratio varies from 1/4 to 4. A selected set of optimization results is displayed in figure 9, and a comprehensive depiction is provided in figures S13-S16. The figures underline the pronounced dependency of optimized designs on the placement of piezoresistors. As the piezoresistor migrates away from the center, the double-cantilever configuration deforms in tandem with increasing piezoresistor aspect ratio, and the degree of deformation is contingent upon the direction of the offset. When the piezoresistor is proximal to the lower or upper edges, intricate new configurations emerge, characterized by holes at the clamped end and within the cantilever. These designs leverage the stress concentration at the corners to maximize sensitivity.
To further investigate the interplay between piezoresistor position and its shape, and their subsequent impact on sensor sensitivity, we conducted a comprehensive evaluation of the two parameters. Our evaluation, as depicted in figure 10(a), indicates that the sensitivity of optimized designs is a nonlinear function of aspect ratio. Moreover, the specific dependency appears to depend stochastically on the piezoresistor position. The optimal aspect ratio of 1/2.5 at zero-offset becomes suboptimal at other piezoresistor positions. The seemingly random patterns in figure 10(a) could be attributed to the solutions being trapped at local minima due to the multidimensional nature and the non-convexity inherent in the current optimization problem. Nevertheless, some general trends are discernable. Specifically, a positive offset typically yields higher sensitivity compared to a negative offset, irrespective of the piezoresistor shape. Furthermore, sensitivity generally decreases with increasing piezoresistor aspect ratios, likely due to the degrading design quality as pointed out in section 3.2. Therefore, a piezoresistor that is short and wide, positioned at the upper corner, will generally lead to optimized designs with superior sensitivity.
For comparative purposes, the sensitivity of rectangular benchmarks with the same piezoresistor configuration is evaluated in figure 10(b). The distribution of piezoresistive signals is shown in figure S17. The sensitivity of these rectangular benchmarks exhibits a distinct dependence on the shape and position of the piezoresistor. Notably, the highest sensitivity is attained at zero offsets, as opposed to the 110 μm offset in figure 10(a). The sensitivity diminishes as the piezoresistor aspect ratio increases, and the decreasing trends exhibit far more regularity compared to the previous case. It is worth noting that an optimal aspect ratio of 1/3 is observed when the offset is −110 and 110 μm, likely due to the optimal alignment of the piezoresistor with the stress profile in bending cantilevers, as can be observed in figure S17. figures S18 and S19 present a detailed breakdown of sensitivity as the weighted sum of stresses in the x and y directions, depending on the position and aspect ratio of the piezoresistor, thereby illustrating the significantly different stress behaviors in optimized designs and rectangular benchmarks.

A simplified model
The inherent complexity of topology optimization makes it challenging to comprehensively analyze the functional consequences of individual geometric components. To investigate the specific contributions from various geometric features-particularly the asymmetric designs observed in the optimized structures-we have developed an FEA model with simplified geometry abstracted from the optimized designs. This model, depicted in figure 11, consists of two rectangular cantilevers joined by a broad base, with the piezoresistor situated at the center. The model is meshed with the same density of quadratic elements as the benchmark models from the previous sections, ensuring comparable accuracy. Utilizing this simplified model, we can systematically manipulate the geometric parameters to examine the origins of high sensitivity and asymmetry resulting from topology optimization. The aspect ratio of the piezoresistor, the lengths of the upper (L upper ) and lower (L lower ) cantilevers, and the gap size (W gap ) serve as control parameters in this investigation.

Effects of gap size
Being a key component in the double-cantilever configuration, we first investigate the influence of the gap between two cantilevers and the dependency of sensitivity on the gap size W .
gap We vary W gap from 0 to 150 μm, maintaining the aspect ratio of the piezoresistor at 1/8 and 1/2.5 for comparison, and setting L upper and L lower at the maximum length of 300 μm. The distributions of piezoresistive signals for a series of geometries are depicted in figure 12, with the stress distributions shown in figures S20 and S21. Without a gap, the stress profile mirrors that of a cantilever under biaxial surface stress, as demonstrated in prior studies [23]. The introduction of a gap, effectively divides one cantilever into two smaller ones, creating a common corner that facilitates a highsensitivity region-an ideal location for the piezoresistor. As the gap expands, the stress concentration lessens, leading to a decrease in sensitivity. Importantly, the shape of the piezoresistor significantly influences sensitivity by determining the extent to which stress concentration generated by the gap is converted into piezoresistive signals. Hence, the relative size of the piezoresistor and the gap may be a critical factor in maximizing sensitivity.
The dependence of sensitivity on gap size varies according to the piezoresistor aspect ratio. In an effort to reproduce the optimized gap sizes in a series of optimized designs showcased in figure 6, we vary both the gap size and the aspect ratio of the piezoresistor to quantitatively investigate their interaction. Figure 13 illustrates a map of sensitivity as a function of gap size and piezoresistor aspect ratio, while figure S22 depicts the corresponding stresses in the x and y directions. In alignment with the previous remarks on figure 8, the sensitivity of the simplified model first increases and then decreases, reaching peak sensitivity at an aspect ratio of 1/2.5 as the piezoresistor aspect ratio increases. Moreover, an optimal gap size exists for all piezoresistor shapes. For lower aspect ratios, the optimal gap size extends up to 80 μm. However, as the aspect ratio increases, the optimal gap size rapidly shrinks and stabilizes around 30 μm, coinciding with the gap sizes measured in optimized designs generated via topology optimization (figure 6), as marked by the black squares in figure 13. It is noteworthy that the dependencies of sensitivity and stresses in different directions in 〈111〉 orientations mirror those in 〈110〉 orientations, albeit with different magnitudes, as demonstrated in figures S23 and S24. Furthermore, to evaluate how accurately the simplified model replicates the sensing behavior of optimized designs obtained from topology optimization, we assess the sensitivity of the simplified model using the same gap sizes as in the optimized designs. These results are presented in figure 14, with the distributions of piezoresistive signals and stresses illustrated in figures S25-S27. Figure 14 highlights the close resemblance between the performance of the simplified model and the optimized designs, barring those that with highaspect-ratio piezoresistors. The discrepancy can be accounted for by the degraded design quality at high piezoresistor aspect ratios, as discussed in section 3.2, which results in increased xx s and , yy s as depicted in figure  S28. Interestingly, the simplified model outperforms the optimized designs at low piezoresistor aspect ratios. This superior performance may be ascribed to the sharp corners in the simplified model, as the minimum feature size in topology optimization is constrained by the dimensions of the mesh and filter. This phenomenon reduces the xx s whilst barely influencing , yy s as indicated in figure S28. Regardless, the strong correspondence between the simplified model and the optimized designs attests to the efficacy of the simplified model in capturing the essential features of optimized designs, thereby enabling further exploration of the impacts of other geometric features.

Effects of asymmetry
The preceding sections have demonstrated that the optimized designs obtained via topology optimization exhibit asymmetrical geometries, with the placement of the piezoresistor significantly influencing the optimized  structures. The asymmetry in optimized designs suggests that the structural materials in the upper and lower halves of the design domain may contribute differently to the sensor performance.
To explore this effect, we maintain the piezoresistor aspect ratio at 1/2.5 and the gap size at 30 μm, and vary the length of one cantilever (e.g., the upper cantilever L upper ) while keeping the other (e.g., L lower ) constant. The resulting distributions of piezoresistive signals and stresses are illustrated in figure 15 and figures S29-S30, respectively. This parametric study is conducted in both 〈111〉 and 〈110〉 orientations, with a quantitative evaluation of sensitivity depicted and compared in figure 16. The sensitivity appears to rise, then plateau, in relation to the length of cantilevers. The saturation sensitivity reaches 13.7 10 4 -(N·m -1 ) -1 , coinciding with the sensitivity of the optimized design at the optimal aspect ratio, as indicated by the dotted horizontal line in the figure. This sensitivity is 23.4% higher than that of the 〈110〉 orientations, which stands at 11.1 10 4 -(N·m -1 ) -1 , thereby demonstrating the superior performance of 〈111〉 orientations even in this optimized geometry.
The saturation length depends on the position of the piezoresistor and aligns well with the average cantilever length measured in the optimized design with a piezoresistor aspect ratio of 1/2.5 (figure 6), as indicated by the vertical lines. The existence of a saturation length elucidates the phenomenon whereby optimized designs with a double-cantilever configuration tend to extend the cantilevers until the sensitivity gain no longer counterbalances the penalty imposed by volume constraints, as discussed in section 3.2.  Another key observation from figure 16 is that when the L upper is fixed and L lower is adjusted, the initial sensitivity is lower, but it reaches a higher saturated sensitivity compared to the inverse case. This feature indicates that the lower cantilever contributes more significantly to sensitivity than the upper one, consistent with the asymmetrical designs obtained thus far and the superior performance of designs that position the piezoresistor at the upper edge of the design domain. This asymmetry is exclusive to 〈111〉 orientations, while it is absent in the conventional 〈110〉 orientations, as evidenced by the completely overlapping curves in figure 16. Furthermore, the asymmetry of sensitivity contribution between the two cantilevers is exacerbated for piezoresistors with high aspect ratios, as illustrated in figure S31.
It should be noted that the preference for the lower cantilever could be trivially switched to the upper cantilever by reversing the third Euler angle used for tensor rotation, specifically, from ( /2, p /4, p

--
). This change does not alter the piezoresistance coefficients but does modify the sign of certain shear components in the stiffness matrix, leading to optimized designs that favor the upper cantilever. The shift of symmetry under tensor rotation allows us to conclude that the asymmetry observed in various optimized designs and the preference for a positive piezoresistor offset are consequences of the anisotropic elasticity of 〈111〉 orientations.

Conclusion
In this study, we leverage topology optimization to explore the optimal design of piezoresistive nanomechanical sensors in highly sensitive 〈111〉 orientations. The efficacy of these designs is rigorously evaluated and benchmarked against traditional rectangular designs, affirming the validity of our approach. Moreover, we develop a simplified FEA model to investigate the functional roles of various geometrical attributes and elucidate the origin of the prevalent asymmetric designs found in optimized designs.
Key findings from our investigation are twofold: Firstly, the optimal design is significantly influenced by the immobilization layer materials. For instance, a double-cantilever design is optimal for thin, soft immobilization layers, whereas a suspended platform design is most suitable for thick, stiff immobilization layers. This transition in optimal design is attributed to the shift in the neutral plane and the dominating stress relaxation mechanism. However, within a technically important parameter range, the double-cantilever design with a low piezoresistor aspect ratio emerges as the most performant design. This design is 22.3% and 23.4% more sensitive than the rectangular benchmarks in 〈111〉 orientations and the same double-cantilever design in 〈110〉 orientations, respectively. This result underscores the superiority of both the design strategy and the material choice.
Secondly, the piezoresistor geometry plays a consequential role in defining the optimized designs and the accompanying sensitivity enhancement. The optimal piezoresistor aspect ratios are found to be asymmetrically position-dependent as a consequence of anisotropic elasticity in 〈111〉 orientations. Hence. the sensitivity Figure 16. The sensitivity of the simplified model under a fixed length of upper/lower cantilever while the length of the other cantilever is varied. The simulation is run twice for 〈110〉 and 〈111〉 orientations. The optimized cantilever lengths measured from topology optimization (aspect ratio 1/2.5 in figure 6) are marked by vertical dashed lines. The sensitivity of the corresponding optimized design is marked by the horizontal dotted line. enhancements obtained in the study are of exemplary value. The actual benefit derived from optimized designs and alternative silicon orientations would depend on the specific application and the fabrication process of the piezoresistive sensor.
Crucially, while the theoretical implications of our findings are profound, it is worth noting the practical challenges associated with their implementation. One notable difficulty in fabricating MEMS devices from 〈111〉 orientations on (110) silicon compared to (100) silicon is the anisotropic etching of silicon. This process, which is dependent on the crystallographic orientation of the silicon wafer, results in different etch rates and shapes for different planes. The {111} planes, for instance, are the slowest etching plane in all anisotropic etchants, and as such, prolonged etching invariably leads to the appearance of {111} facets at the sidewalls of the fabricated structures [60]. This occurrence can influence the performance of MEMS devices that rely on accurate geometries and dimensions.
Despite the potential challenges with fabrication, this study provides invaluable guidelines for future design considerations utilizing alternative silicon orientations. This work can potentially improve the capabilities of piezoresistive nanomechanical sensors across a wide array of disciplines, including but not limited to healthcare, agriculture, and homeland security.