Investigation of polyelectrolyte gel based electronic devicesbyVasilii TriandafilidiBSc, Moscow Institute of Physics and Technology, 2013MSc, University of British Columbia, 2015A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Chemical and Biological Engineering)The University of British Columbia(Vancouver)April 2020c© Vasilii Triandafilidi, 2020The following individuals certify that they have read, and recommend to the Fac-ulty of Graduate and Postdoctoral Studies for acceptance, the thesis entitled:Investigation of polyelectrolyte gel based electronic devicessubmitted by Vasilii Triandafilidi in partial fulfillment of the requirements for thedegree of Doctor of Philosophy in Chemical and Biological Engineering.Examining Committee:Dr. Savvas G. Hatzikiriakos, Chemical and Biological EngineeringSupervisorDr. Joerg Rottler, Physics and AstronomyCo-supervisorDr. James Feng, Chemical and Biological EngineeringSupervisory Committee MemberDr. John Madden, Electrical and Computer EngineeringSupervisory Committee MemberDr. Stavros Avramidis, Faculty of ForestryUniversity ExaminerDr. Peter Englezos, Department of Chemical and Biological EngineeringUniversity ExaminerDr. Vlasis G. Mavrantzas, Department of Chemical Engineering Laboratory ofStatistical Thermodynamics and Macromolecules at ETH ZurichExternal ExamineriiAbstractWe present a molecular dynamics (MD) and experimental study of polyelectrolyte(PE) gel based electronic devices such as sensors and diodes. We first perform anMD study of two PE gels with different degrees of ionization coupled in a slabgeometry. Our simulations show that a pressure gradient emerges between the twogels that results in the buildup of a Nernst-Donnan potential. The Nernst-Donnanpotential at the interface is found to scale linearly with temperature with the coeffi-cient of proportionality given by the fraction of concentrations of the uncondensedcounterions. We show that the potential difference can also be expressed as a lin-ear function of the lateral pressure, thus providing a molecular interpretation ofthe piezo-ionic effect. These findings provide further insight onto the behaviour ofsoft-sensors in the equilibrium regime with no salt ion/solvent fluxes.We also perform an MD study of a junction of two oppositely charged PEnetworks, and compare the ion densities and electrostatic field to a correspondingcontinuum Poisson-Boltzmann (PB) model. At low electrostatic coupling strength,the PB model reproduces the MD simulation results for density and electric fieldthroughout the gel very well. At higher electrostatic coupling and higher degreesof ionization, the standard PB fails to predict the MD profiles at the diode inter-face due to counterion condensation, network collapse and field-induced gel de-formation. In fact, MD simulations predict that the rectifying behavior of diodesoperating in such regimes will be much reduced. We develop a modified PB modelthat accounts for these effects, show that it produces better agreement with theMD results, and can be used for improved modeling of Polyelectrolyte Gel Diodes(PGDs).Additionally, we perform a systematic stress-test of the predictions of the Ya-iiimamoto and Doi [77] theory of the PGDs under linear sweep and step bias. Wehave found that the predictions of the functional forms of current-voltage (I-V)curves of the Yamamoto theory hold well. They predict an exponential increase inthe regime of the forward bias as well as the square-root dependency for the regimeof the reverse voltage.ivLay SummaryWith a rapid advancement of robotics, many researchers are becoming increas-ingly interested in applying electronic devices in the medical field, producing bio-inspired mechanisms. However, this requires materials that are bio-compatible andhave good mechanical properties like flexibility and stretchability. Unfortunately,most materials used in modern day electronics are either solid like in electric wiresor liquid like in car electrolytes. Polyelectrolyte (PE) hydrogels are a promis-ing material that combine the electrical properties of the incumbent materials withsuperior mechanical properties. Despite that they remain underused due to theirnovelty and our lack of understanding of their structure. The current work focuseson systematic investigation of PE hydrogels using computer simulation techniquesas well as experiments. We focus on studying two most fundamental electronicdevices sensors and diodes, and study their analogues comprising of PE hydrogels.Our work provides more insights into the nanoscopic structure of hydrogels andallow for better design of the soft electronic devices.vPrefaceThe material in this thesis is the original work of the author, Vasilii Triandafilidi(VT). The investigation described in Chapter 4 and Chapter 5 including identifi-cation of the problem, design of the computational experiments, computer simula-tions and analysis of the data is conducted solely by VT. Several portions of thisthesis have been published as original manuscripts in peer-reviewed journals withVasilii Triandafilidi as the primary author: Chapter 4 is published in Triandafilidiet al., Soft Matter (2018), a version of Chapter 5 in Triandafilidi et al., Soft Matter(2020). Here problem identification, simulation, and analysis of the data is con-ducted solely by VT; manuscripts were edited by Dr. Joerg Rottler and Dr. SavvasG. Hatzikiriakos.Experimental observations reported Chapter 6 are conducted in conjunctionwith two other University of British Columbia graduate students Parya Keyvani(PK), Kudzi Nayamayaro (KN). The synthesis of the structures was done by PKand KN, the electrochemical measurements were done by PK and VT. The datawas compiled, analyzed and matched with theoretical models by VT. A versionof Chapter 6 is currently being submitted for publication. All other chapters areoriginal work, first published in this document.All work was supervised by Dr. Savvas G. Hatzikiriakos and Dr. Joerg Rot-tler including guidance in analyses and project direction. Resources of ComputeCanada have been used to run the simulations.viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1.1 New materials for new applications . . . . . . . . . . . . 11.1.2 Polyelectrolyte hydrogels as soft sensors . . . . . . . . . 31.1.3 Polyelectrolyte Gel Diodes . . . . . . . . . . . . . . . . . 51.2 Piezo-ionic sensors . . . . . . . . . . . . . . . . . . . . . . . . . 51.2.1 Literature review . . . . . . . . . . . . . . . . . . . . . . 51.2.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . 121.2.3 Outline of the work on the piezo-ionic effect . . . . . . . 121.3 Polyelectrolyte Gel Diodes (PGDs) . . . . . . . . . . . . . . . . 13vii1.3.1 Literature review . . . . . . . . . . . . . . . . . . . . . . 131.3.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . 161.3.3 Outline of PGD work . . . . . . . . . . . . . . . . . . . . 171.4 General outline of the thesis . . . . . . . . . . . . . . . . . . . . 182 Theoretical background . . . . . . . . . . . . . . . . . . . . . . . . . 202.1 Theoretical aspects of polyelectrolyte gels . . . . . . . . . . . . . 202.1.1 A molecular view of polyelectrolytes . . . . . . . . . . . 202.1.2 Introduction to Polymer Physics . . . . . . . . . . . . . . 222.1.3 Counterion condensation . . . . . . . . . . . . . . . . . . 282.1.4 Donnan equilibrium . . . . . . . . . . . . . . . . . . . . 292.2 Theory of Polyelectrolyte Gel Diodes . . . . . . . . . . . . . . . 302.2.1 Poisson Boltzmann equation . . . . . . . . . . . . . . . . 302.2.2 Transport equations . . . . . . . . . . . . . . . . . . . . . 312.2.3 Yamamoto and Doi theory of PE diodes . . . . . . . . . . 323 Methodology and materials . . . . . . . . . . . . . . . . . . . . . . . 403.1 Method of Molecular Dynamics (MD) . . . . . . . . . . . . . . . 403.1.1 Introduction: What is MD? . . . . . . . . . . . . . . . . . 403.1.2 Core of MD modeling: equations of motion . . . . . . . . 423.1.3 The ”nuts and bolts” of the MD simulation . . . . . . . . 433.1.4 Details of implementation . . . . . . . . . . . . . . . . . 523.2 Continuum theory and modelling . . . . . . . . . . . . . . . . . . 563.3 Experimental methods . . . . . . . . . . . . . . . . . . . . . . . 573.3.1 Fabrication of PGD . . . . . . . . . . . . . . . . . . . . . 574 Piezoionic effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 624.1 Three regimes of Donnan equilibrium . . . . . . . . . . . . . . . 624.2 Potential energies . . . . . . . . . . . . . . . . . . . . . . . . . . 664.3 Nernst-Donnan model for the electrostatic potential . . . . . . . . 674.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 705 Molecular Dynamics modeling of PE diodes . . . . . . . . . . . . . . 715.1 Weakly ionized chain with f = 0.1 . . . . . . . . . . . . . . . . . 71viii5.1.1 Weak coupling (lb = 1) . . . . . . . . . . . . . . . . . . . 715.1.2 Strong coupling (lb = 5) . . . . . . . . . . . . . . . . . . 775.2 Effect of ionization ( f = 0.25) . . . . . . . . . . . . . . . . . . . 805.3 Modified Poisson Boltzmann model . . . . . . . . . . . . . . . . 835.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 856 Experimental testing of Yamamoto and Doi theory . . . . . . . . . . 876.1 Investigation of the J(V) curves under linear sweep of voltages . . 876.1.1 Voltammograms under linear sweep . . . . . . . . . . . . 876.1.2 Regime of forward bias . . . . . . . . . . . . . . . . . . . 886.1.3 Regime of backward bias . . . . . . . . . . . . . . . . . . 906.2 PGDs under step voltage . . . . . . . . . . . . . . . . . . . . . . 916.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 927 Conclusion and future work . . . . . . . . . . . . . . . . . . . . . . 947.1 Piezo-ionic effect . . . . . . . . . . . . . . . . . . . . . . . . . . 947.2 MD simulation of PE diode . . . . . . . . . . . . . . . . . . . . . 957.3 Experimental investigation of PE diode . . . . . . . . . . . . . . 977.4 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99A Supporting Materials . . . . . . . . . . . . . . . . . . . . . . . . . . 108A.1 Validation of the MD model . . . . . . . . . . . . . . . . . . . . 108ixList of TablesTable 3.1 Table of conversion between real and LJ units . . . . . . . . . 47Table 3.2 A summary of the simulation parameters used in the present work. 55Table 3.3 Zeta potential, crystallinity index, charge density from conduc-tometric titration and charge density from elemental analysisof unmodified CNC, nCNC and pCNC. Error bars represent astandart deviation from the mean after 3 measurements . . . . 59xList of FiguresFigure 1.1 Schematic representation of a PE hydrogel . . . . . . . . . . . 4Figure 1.2 A schematic representation of the explanation of the piezo-ionic effect proposed by Sawahata et al. [50, 64] (see text). . 7Figure 1.3 A schematic representation of the explanation of the piezo-ionic effect proposed by Sawahata et al. [50, 64] (see text). . 8Figure 1.4 A: A molecular picture of the Nernst-Donnan equilibrium be-tween two gels with different ionizations studied in experi-ments. Gels with a higher degree of negative backbone ion-ization have an intrinsically lower negative electrical poten-tial when compared to the gels with lower degree of negativebackbone ionization. The system consists of negative back-bone ions (white), floating counterions (green), dissolved saltions (orange and purple) and water molecules. B: A schematicrepresentation of the explanation of the piezo-ionic effect pro-posed by Sawahata et al. [50, 64] (see text). . . . . . . . . . . 10xiFigure 1.5 A: Schematic representation of the interplay between the os-motic pressure Πosm and elastic energy of the network Πel inthe Donnan equilibrium. B: Schematic illustration of coun-terion condensation. When the Bjerrum length lB of the sys-tem becomes bigger than the distance between two neighbor-ing backbone macroions a≤ lB, the counterions become boundto the backbone ions and stop contributing to the osmotic pres-sure. The system consists of negative backbone ions (white),floating counterions (green). . . . . . . . . . . . . . . . . . . 11Figure 1.6 Schematic representation of the polyelectrolyte gel diode in Arectifying and B blocking mode. Due to applied voltage, wa-ter present in the gel diode undergoes an electrolysis processgenerating H+ and OH– ions at the cathode and anode, respec-tively and, simultaneously, produces electrons in the externalcircuit. The asymmetry of the potential drop either attracts thecounterions to the interface (thus allowing OH– , H+ groupsto recombine and rectify the current) or pushes counterions toelectrodes (thus blocking OH– , H+ groups from recombinationand blocking the current). . . . . . . . . . . . . . . . . . . . 13Figure 2.1 An example of a flexible PE. Left: constitution formula forsulfonated polystyrene and its sodium counterions; right: acoarse-grained representation that forsakes the chemical detailto focus only on the physics of the system (taken from [45]) . 21Figure 2.2 An example of a rod-like PE. Left: constitution formula forpoly(paraphenylene) with iodine counterions; right: a coarse-grained representation (taken from [45]) . . . . . . . . . . . . 22Figure 2.3 A schematic representation of an ideal chain (taken from [59]) 23Figure 2.4 Left: a schematic representation of thermal blobs in a polymerunder tension; middle: chain extensions for different modelsof polymer chains; right: comparison of a chain extension be-tween WLC model and experimental results for λ -DNA (takenfrom [59]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24xiiFigure 2.5 A schematic representation of an ideal chain (top) vs real chain(bottom) of the same counter length Rmax = Nmb stretched bythe same force−→fx (taken from [59]) . . . . . . . . . . . . . . 26Figure 2.6 Electrostatic potentials in the equilibrium. Polyelectrolyte geldiodes (PGD) are the double layers of two oppositely chargedpolyelectrolyte gels that are sandwiched by two symmetric elec-trodes. Experiments have shown that PGDs show very smallelectric currents for reverse applied voltages (a), V0 < 0, andshow very large electric currents for forward applied voltages(b), V0 > 0, for the cases in which salt impurity is eliminated[10, 82]. The bottom figures schematically show electrostaticpotentials ϕ(z) (the black solid curve) and the local densitiesof positive and negative counterions, c+(z) and c−(z), (the redand blue solid curves) in a half of the PGD, z > 0, for the casesin which electrochemical reactions at the electrodes are sup-pressed. The ionic groups of polyelectrolyte gels are treated asuniform fixed electric charges (these charges are shown by +and - with filled circles) and the distributions of their counteri-ons are functions of electrostatic potentials (these counterionsare shown by + and - with open circles). PGDs show deple-tion layers, where counterions are depleted, at the interfaceregion between the two oppositely charged gels, 0 < z < zd ,and the surfaces of the electrodes, ze < z < L/2. Electrostaticpotentials ϕ(z) are approximately constant cf in the bulk re-gion of polyelectrolyte gels. Bulk electrostatic potentials cfare (halves of) potential drops across the depletion layer atthe interface between the two gels and thus are generated bythe fixed charges of the polyelectrolyte gel in this layer (takenfrom Yamamoto and Doi [77]) . . . . . . . . . . . . . . . . . 38xiiiFigure 2.7 Voltammogram (I(V)) profiles of PGD in the system withoutadditional salt. Electric current densities Jele =(= eJH =−eJOH)(that are rescaled by exchange current density eJ0, see Meth-ods) are calculated as functions of applied voltages V0 (in theunit of the thermal voltage) for reverse (a) and forward volt-ages (b). The number of hydrogen and hydroxide ions that areaccumulated in the bulk of the gel increases with decreasingthe values of k f 0 (that is defined by equation (2)). We usedk f 0 = 1.0 ·106 (black solid curve), 5.0 ·106 (blue) and 1.0 ·107(red) for the calculations. The value of (rescaled) gel thick-ness k f L is fixed to 1.0 · 104. The equations that are used forthe calculations are shown in Section 2.2.3, 2.2.3. We useda simple scaling argument (shown in Methods) to predict the(black) broken curve in a. The (black) broken curve in b showsthe cases in which all of the applied voltages are distributed tothe surfaces of the electrodes (taken from Yamamoto and Doi[77]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39Figure 3.1 The length and time scales of MD simulations compared toavailable experimental techniques and other theoretical meth-ods (courtesy of Kochmann research group [32]). . . . . . . . 41Figure 3.2 The algorithm of an MD simulation. First, initial configura-tions are created for every atom in the system. Then forcesapplied to every atom are calculated by setting the interaction”strength” through inter-atomic potential. Then equations ofmotion (Equation 3.1) are integrated. . . . . . . . . . . . . . 43Figure 3.3 The coarse-graining methodology where the whole C−H3 unitsare represented as a single bead. . . . . . . . . . . . . . . . . 46Figure 3.4 The potential profile for MD simulation. . . . . . . . . . . . . 48Figure 3.5 Schematic representation of the periodic boundary conditions,where the systems boundaries are shrink wrapped to avoid thefinite size effects. . . . . . . . . . . . . . . . . . . . . . . . . 50Figure 3.6 Schematic representation of the p-type gel synthesis. . . . . . 57xivFigure 3.7 Schematic representation of the gel synthesis. . . . . . . . . . 57Figure 3.8 Schematic representation of the electrode setup. . . . . . . . . 58Figure 3.9 Zeta potential from conductometric titration of unmodified CNC,nCNC and pCNC. Error bars represent a standard deviationfrom the mean after 3 measurements. . . . . . . . . . . . . . 60Figure 3.10 Schematic representation of the fabrication method to fabricatea PGD. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60Figure 3.11 Schematic representation of study protocol of a PGD. A lin-ear sweep between voltage range at different rate(black), and astep bias(red). . . . . . . . . . . . . . . . . . . . . . . . . . 61Figure 4.1 Snapshots of a system with ionization degrees f1 = 0.1 andf2 = 1.0 at three temperatures T = 1.1,0.3,0.1 (at lB = 0.9,3.3,10).In snapshots white balls are: uncharged backbone monomers,red: negatively charged monomers, blue: tethered static monomers,green: positively charged counterions. For illustration pur-poses, only counterions within the dashed box are shown. Regime1 (weak-weak): weak coupling for both sides of the sensorwith no counterion condensation on both sides Regime 2 (weak-strong): weak coupling on the left side, and strong coupling onthe right side of the sensor Regime 3 (strong-strong): strongcoupling regime for both parts of the system . . . . . . . . . 63Figure 4.2 Analysis of the fraction of free counterions (open symbols)compared to the Manning theory (filled symbols) for the gelslabs with ionization degrees f1 = 0.1 and f2 = 1.0 (blue) andf1 = 0.25 and f2 = 0.5 (cyan). A: Fraction of free counterionsf f ree,1 = f1(1−Q1) on the left side of the gel. B: Fraction offree counterions f f ree,2 = f2(1−Q2) on the right side of thegel. C: Ratio f f ree,2/ f f ree,1. Vertical dashed lines correspondto the temperatures where counterion condensation occurs inone part of the gel as indicated by the change in the curve be-haviour. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65xvFigure 4.3 Normalized energy profiles (A: Coulombic contribution, B:Lennard-Jones contribution) of the system with f1 = 0.1, f2 =1.0 for different regimes. Here N is the total number of atomsin the system. . . . . . . . . . . . . . . . . . . . . . . . . . 67Figure 4.4 Counterion concentration profiles of the system with f1 = 0.1,f2 = 1.0 for different temperature regimes. Colors are the sameas in Figure 4.3 . . . . . . . . . . . . . . . . . . . . . . . . . 68Figure 4.5 Normalized electrostatic potential as a function of temperature.Circles: q∆ϕ∗ = q∆ϕ/ ln(cfree,2/cfree,1). Red: f1 = 0.1, f2 =0.2; magenta: f1 = 0.1, f2 = 0.3; green: f1 = 0.1, f2 = 0.5;blue: f1 = 0.1, f2 = 1.0. The dashed line has slope unity. . . 69Figure 4.6 Potential difference as a function of the lateral pressure differ-ence. Red: f1 = 0.1, f2 = 0.2; magenta: f1 = 0.1, f2 = 0.3;green: f1 = 0.1, f2 = 0.5; blue: f1 = 0.1, f2 = 1.0. Here V isthe volume of the simulation box. . . . . . . . . . . . . . . 70Figure 5.1 Ion spatial density profile of the system with f = 0.1, lb = 1at different strengths of applied electric field E = 0 (top), E =−0.1 (middle), E = 0.1 (bottom). Blue: negatively chargeditinerant ions; red: positively charged itinerant ions; green:gel backbone ions (both positively and negatively charged);solid line: PB model solution; dashed line: MD profiles. Errorbars represent±1σstd confidence intervals (for clarity purposesonly every 6th error bar is shown). . . . . . . . . . . . . . . 72Figure 5.2 Spatial profile of the electric field inside of the diode at differ-ent strengths of applied field (black: E = 0, blue: E = −0.1,red: E = +0.1) for the weakly ionized system in the limit ofweak coupling ( f = 0.1, lb = 1). Dashed lines represent theMD data; solid lines are the corresponding PB model curves.Error bars represent ±1σstd confidence intervals (for claritypurposes only every 6th error bar is shown). . . . . . . . . . . 73xviFigure 5.3 Ion spatial density profile of the system with f = 0.25, lb = 1at different strengths of applied electric field E = 0 (top), E =−0.1 (middle), E = 0.1 (bottom). Blue: negatively chargeditinerant ions; red: positively charged itinerant ions; green:gel backbone ions (both positively and negatively charged);solid line: PB model solution; dashed line: MD profiles. Errorbars represent±1σstd confidence intervals (for clarity purposesonly every 6th error bar is shown). . . . . . . . . . . . . . . . 75Figure 5.4 Ion spatial density profile of the system with f = 0.1, lb = 5at different strengths of applied electric field E = 0 (top), E =−0.1 (middle), E = 0.1 (bottom). Blue: negatively chargeditinerant ions; red: positively charged itinerant ions; green:gel backbone ions (both positively and negatively charged);solid line: PB model solution; dashed line: MD profiles. Errorbars represent±1σstd confidence intervals (for clarity purposesonly every 6th error bar is shown). . . . . . . . . . . . . . . . 79Figure 5.5 Spatial profile of the electric field inside of the diode at differ-ent strengths of applied field (black: E = 0, blue: E = −0.1,red: E = +0.1) of the weakly ionized system in the limit ofstrong coupling ( f = 0.1, lb = 5). Dashed lines represent theMD data; solid lines are the corresponding PB model curves.Error bars represent ±1σstd confidence intervals (for claritypurposes only every 6th error bar is shown). . . . . . . . . . . 80Figure 5.6 Ion spatial density profile of the system with f = 0.25, lb = 1at different strengths of applied electric field E = 0 (top), E =−0.1 (middle), E = 0.1 (bottom). Blue: negatively chargeditinerant ions; red: positively charged itinerant ions; green:gel backbone ions (both positively and negatively charged);solid line: PB model solution; dashed line: MD profiles. Errorbars represent±1σstd confidence intervals (for clarity purposesonly every 6th error bar is shown). . . . . . . . . . . . . . . . 81xviiFigure 5.7 Spatial profile of the electric field inside of the diode at differ-ent strengths of applied field (black: E = 0, blue: E = −0.1,red: E = +0.1) in the limit of weak coupling (lb = 1) andionization fraction f = 0.25. Dashed lines represent the MDdata; solid lines are the corresponding PB model curves. Errorbars represent±1σstd confidence intervals (for clarity purposesonly every 6th error bar is shown). . . . . . . . . . . . . . . . 82Figure 5.8 Ion spatial density profile of the system with f = 0.25, lb = 5at different strengths of applied electric field E = 0 (top), E =−0.1 (middle), E = 0.1 (bottom). Blue: negatively chargeditinerant ions; red: positively charged itinerant ions; green:gel backbone ions (both positively and negatively charged);solid line: PB model solution; dashed line: MD profiles. Errorbars represent±1σstd confidence intervals (for clarity purposesonly every 6th error bar is shown). . . . . . . . . . . . . . . . 83Figure 5.9 Spatial profile of the electric field inside of the diode at differ-ent strengths of applied field (black: E = 0, blue: E = −0.1,red: E = +0.1) in the limit of weak coupling (lb = 5) andionization fraction f = 0.25. Dashed lines represent the MDdata; solid lines are the corresponding PB model curves. Errorbars represent±1σstd confidence intervals (for clarity purposesonly every 6th error bar is shown). . . . . . . . . . . . . . . . 84Figure 5.10 Spatial profile of the electric field inside of the diode at differ-ent ionizations ( f = 0.1 top, f = 0.25 bottom) and couplings(black: lb = 5, red: lb = 1). Circles represent the MD data;solid lines are the corresponding modified PB model curveswith parameter d described in Table 3.2. Error bars represent±1σstd confidence intervals (for clarity purposes only every6th error bar is shown). . . . . . . . . . . . . . . . . . . . . . 85Figure 6.1 The resulting IV curves obtained by sweeping through voltages-3 to 3 [V] with the rate of 50 [mV/s]. . . . . . . . . . . . . . 88xviiiFigure 6.2 Lines: the logarithm of the electric current density and the cor-responding data fit for the regime of the forward bias. Thefindings concur with the theoretical prediction of Yamamotoand Doi [77] with J ∼ exp(−∆V ). Results are obtained bysweeping through voltages -3 to 3 [V] with different sweeprates). Circles: the resulting IV curves obtained by taking ter-minal values of currents in forward and backward bias regimesunder different magnitudes of applied voltage . . . . . . . . . 89Figure 6.3 Electric current density and the corresponding data fit for theregime of the backward bias. The findings concur with the the-oretical prediction of Yamamoto and Doi [77] with J ∼ V 1/2.Results are obtained by sweeping through voltages -3 to 3 [V]with the rate of 50 [mV/s]). . . . . . . . . . . . . . . . . . . . 90Figure 6.4 The response of a PGD to a step-wise voltage input. Green linesindicate the terminal value of the currents (in the backward andforward regimes) . . . . . . . . . . . . . . . . . . . . . . . . 91Figure A.1 An electrolyte system placed in a box with reflective sides andperiodic boundary conditions in lateral directions and corre-sponding electrical field profiles. . . . . . . . . . . . . . . . . 109xixGlossaryFENE Finite-Extensible-Nonlinear-ElasticLJ Lennard-JonesNPT A isothermal–isobaric ensemble where the total number of particles inthe system (symbol: N), the system’s pressure (symbol: P), as well as thetemperature of the system is preserved (symbol: T)NVE A microcanonical ensemble where the total number of particles in thesystem (symbol: N), the system’s volume (symbol: V), as well as thetotal energy in the system (symbol: E)NVT A canonical ensemble where the total number of particles in the system(symbol: N), the system’s volume (symbol: V), as well as thetemperature of the system is preserved (symbol: T)PAA Poly (Acrylic Acid)PAM Poly (Arcylamide)PB Poisson BoltzmannPBC Periodic boundary conditionsPE PolyelectrolytePGD Polyelectrolyte Gel DiodePMAA Poly (Methacrylic Acid)WLC Worm like chain model of a polymerxxAcknowledgmentsI thank my fellow labmates 1 Tanja Tomkovic, Vinod Konaganti, Marzieh Ebrahimiand many others for the stimulating discussions, for the sleepless nights we wereworking together before deadlines, and for all the fun we have had in the last fouryears.The initial inspiration for this research was born from tireless conversationswith Dr. Saquib Sarwar and the Mechatronics group2, Claire Preston, Yuta Dobashi,Dr. John Madden and others. They provided me an opportunity to join their team,and gave access to the laboratory and research facilities. Without their precioussupport it would not be possible to conduct this research. Many thanks to the SoftMatter group 3 for insightful conversations on Molecular Dynamics and our semi-nars where we combined science with tasting of wine and cheese.I greatly appreciate the help of my committee members Dr. James Feng, Dr.John Madden for their valuable feedback throughout my PhD. I would also like toacknowledge Dr. Felix Blyakmann, Dr. Alexey Dobrynin, Dr. Aykut Erbas, Dr.Tetsuya Yamamoto, Dr. Giovanni Natale for valuable discussions and feedback.Their insightful comments and encouragement, as well as hard questions motivatedme to widen my research horizons from various perspectives.I would like to thank people who helped me start my career in science. MyPhysics and Maths teachers at 31st lyceum, my teachers at MIPT Moscow Inst ofPhys and Technology. Dr. Andrey Gavrikov for helping me through those difficult4 years, Dr. Gleb Sklizkov for providing me with materials on Physics that I use to1Polymer Rheology group, Dept. of Chem Bio Engineering at the University of British Columbia2Molecular Mechatronics group at the University of British Columbia3Soft Condensed Matter And Computational Materials Physics at University of British ColumbiaDept. of Physics and Astronomy and Quantum Matter Institutexxithis day! My undergraduate supervisors Dr. Vladimir Stegailov, Genry Norman atJIHT 4.I express my sincere gratitude to my advisors Prof. Savvas G. Hatzikiriakosand Prof. Joerg Rottler for the continuous support of my Ph.D. study, for their mo-tivation, and immense knowledge. Their patience through reading the early draftsof my papers which now make even me embarrassed. I could not have imaginedhaving better advisors and mentors for my Ph.D. study.Last but not the least, I would like to thank my family. My father for his wis-dom and great conversations that we have had over the last couple of years. Mymother for believing in me when I was just starting to connect my life with sci-ence. My brother and our tireless conversations about cars and electric engines.Big thanks to Galanopoulos, Dordzhiev and Kritharis families for being my fam-ily away from home. For my koumparos Nicholaos Galanopoulos and my friendAthanasios Kritharis, it would be very gray 4 years without you physically andspiritually being near!This research was supported in part by a University of British Columbia FourYear Fellowship, in part by Compute Canada (www.computecanada.ca) and in partby Natural Sciences and Engineering Research Council (NSERC).4Molecular Modeling group at Joint Inst of High Temperatures of Russian Academy of SciencexxiiChapter 1IntroductionAlways start with ”why?” — Simon Sinek1.1 Introduction1.1.1 New materials for new applicationsThe remarkable progress of electronics in the past 50 years has now had its foot-print in almost every modern day industry. What started with a seminal discoveryof a solid-state transistor in 1947, led to creation of a personal computers, musicplayers and mobile phones. Nowadays, we can see electronic devices to be usedeverywhere ranging from robotics, wearables to medical implants [5, 8, 22, 31,35, 50, 62, 72]. Surprisingly though, the growing success of solid state electronicsproduced a demand for new materials that is far greater than what the materials oftoday can offer.With a rapid advancement of robotics, many researchers are becoming increas-ingly interested in applying electronic devices in the medical field, producing bio-inspired or bio-hybrid mechanisms. However, to apply robotics at the interface ofbiology and engineering one needs new materials with advanced mechanical andelectric properties. The requirements are particularly strict as the materials needto be self-powered, bio-compatible and have good mechanical properties like highstretchability and flexibility [26, 28, 34].1Two applications of a particular interest are soft sensors and soft diodes. Sen-sors are devices that are able to change their properties upon application of externalstimuli. An example of a such device would be a piezo-electric sensor that convertsmechanical pressure to an electric signal like in a door buzzer. Other solid sensors(known as capacitive) are widely used, for example in our smartphones or smartwatches. Using those sensors to replace human skin, however, is problematic. Thesolid sensors don’t have the right mechanical properties, are not flexible, often notbio-compatible and can’t feel all types of mechanical impact (such as shear) [62].Thus, it is essential to recreate a material that can improve on the sensing capabili-ties of a ”solid-state” materials with the mechanical properties of a soft human-skinlike matter.Another, particularly interesting application is creating a soft diode. A diodeis a device that can rectify an electrical current. When a positive external voltage(forward bias) is applied it allows the electrical current to pass through it. How-ever, when a negative voltage is applied (reverse bias) it acts as an insulator. Thesolid-state diodes became the cornerstone upon which the solid-state transistorswere built. The transistors in their turn were the building bricks of computers,performing all the necessary logical operations. Decades of innovation allowed usto perfect the creation of solid-state diodes and make them ubiquitous. We havelearned to make devices from these materials to be small and have good electricalproperties like high power output that were necessary for electronic devices. How-ever, they are unable to satisfy the demands for superior mechanical properties suchas elasticity and bio-compatibility.To match the gap between the mechanical and electrical properties of materialspolyelectrolyte (PE) gels have been introduced. Polyelectrolyte gels are more flex-ible than solid materials (can be stretched up to 3-4x their length) and can be easilycontained. They can be synthesized to be non-toxic, bio-friendly and cheap. Be-low we will show the unique properties of these remarkable materials and illustratehow they could be used to bridge the gap between the future demands of electronicindustry and current scope of available materials.21.1.2 Polyelectrolyte hydrogels as soft sensorsGels in a dry form have a neutral backbone and have no ions. They consist ofa backbone-forming polymer, water and dissolved salt molecules (if added). Forsome gels, in the presence of a polar solvent like water, backbone monomers disso-ciate into charged ionized backbone and oppositely charged itinerant (also knownas floating) counterions[1, 50, 59]. These gels are referred to as polyelectrolyte(PE) gels and consist of a backbone like carboxylic group in Poly Acrylic Acid(PAA), the gels with a non-dissociating backbone like Poly Acrylamide (PAM) arecalled neutral.The backbone chains are connected via chemical covalent cross-links, entan-glements or via physical bonds (e.g., hydrogen-bonds) and form a three-dimensionalnetwork that traps any floating ions inside of it. These ions could be the counteri-ons or any additional ions that are introduced into the network like, for example,from dissociation of salt NaCl −−→ Na+ +Cl– . The floating ions allow the gelto transmit an electric current. Just like the electrons are transmitting the currentin metal, counterions (in PE gels) or ions of brine (in neutral gels with dissolvedsalt) transmit the electric signal within the gel. Unlike neutral gels (like PAM),PE gels do not require any additional salt to transmit or generate an electrical sig-nal, which becomes important when different mechanisms sensing are discussed inSection 1.2.The unique properties of hydrogels come from their peculiar structure of asemi-fixed stretachable network as well as a highly hydrophilic backbone produc-ing an osmotic pressure. This pressure absorbs the solvent molecules to solvate thebackbone and swells the network. The swelling continues until the free energy ofnetwork stretching becomes equal to the free energy of the ion solvation. The os-motic pressure provides hydrogels with an enormous amount of applications. Forexample, when immersed in a solution hydrogels can absorb enormous amountsof solvent (100x) their dry size, swelling to large volumes. This allowed hydro-gels to be used in our everyday lives for hygiene products, for drug delivery andmicrofluidics [54], and for water desalination [33].The hydrogels are extremely inexpensive, flexible and can be synthesized torespond to external stimuli like temperature, pH or mechanical stimuli. These3Figure 1.1: Schematic representation of a PE hydrogelstimuli-responsive hydrogels are usually referred to as “smart” or ”intelligent” hy-drogels [12, 61, 64]. Polyelectrolyte hydrogels, also have amazing mechanicalproperties. One could stretch them 3x their size, they have great ductility andcould be synthesized to be bio-compatible. These qualities make PE hydrogels agreat candidate to be used as soft sensors. 11There are of course number of mechanical and electrical disadvantages when using the sensors,which have to be taken into consideration. Hydrogels are materials made of water, hence are subjectto water evaporation and change in its properties. Also, hydrogels have poorer signal propertiesdue to high noise to signal ratio. Besides, the frequency of the signal at which gels can operate atare limited by the ion diffusion and remain in the kHz range. Notwithstanding the disadvantages,elastomer matrices could be used for better evaporation properties, as well as better digital filterscould be used to allow for a better signal-to-noise ratio.41.1.3 Polyelectrolyte Gel DiodesDepending on the sign of the charge of the backbone, polyelectrolyte gels can beeither positive (comprising of positively charged basic monomer units and nega-tively charged counter ions) or negative (comprising of negatively charged acidicmonomer units and positive counter ions) 2.When the gel is sandwiched between electrodes and a voltage is applied, waterpresent in the gel undergoes an electrolysis process generating H+ and OH– ions atthe anode and cathode, respectively. Simultaneously, electrons are generated andtransported in the external circuitry of the diode. The asymmetry of the gel and thepotential drops at the interface give rise to the rectifying ability of the gel 3, i.e.,depending on the sign of the applied voltage, the diode will either block the electriccurrent or rectify it, shown in Figure 1.6.Despite having superior mechanical properties, the PE hydrogels are still notwidely spread. Partially, it is due to their relative novelty, as well as our lack ofunderstanding of physics that drives the sensing and rectifying behavior. The lackof understanding prevents us from being able to fine-tune the materials to reachsuperior performance. Below we will dive deeper into the history of PE hydrogels,and the available literature on using PE hydrogels as soft sensors and diodes aswell as the gaps in our current understanding of them.1.2 Piezo-ionic sensors1.2.1 Literature reviewAs was mentioned above, by choosing the monomers of a polyelectrolyte hydrogel,one can program it to respond to external stimuli like temperature, pH or mechan-ical deformation (as in touch sensors). The appearance of an electrical voltagebetween different regions due to a non-uniform mechanical stress is called piezo-ionic or mechano-electric effect [12, 61, 64]. To optimize the performance of hy-drogels in touch sensors, one needs to understand in detail the motion of ions when2For example, one could use polycationic Polydiallyldimethylammonium chloride (PDADMAC)with negatively charged Cl– counterions and polyanionic polystyrene sulfonate (PSS) with positivelycharged H+ counterions on the other side of the of diode [10]3A more thorough introduction of PGDs is described in Section 2.2.35pressure is applied.From the empirical evidence it seems like the piezo-ionic effect could be aproduct of two different physical phenomena that could be explained via eitherequilibrium Donnan type analysis [23, 25, 54] or non-equilibrium thermodynam-ics [7, 20, 61, 64, 69]. The non-equilibrium regime appears due to salt ion orsolvent fluxes inside of the gel. At high concentrations of dissolved salt C ≈ 1 M,i.e., low Debye lengths (λD = 0.3nm) all electrostatic interactions become screenedand solvents play little role in determining the system response. Hence ion motionbecomes a dominant factor in the gel behaviour.4 In this regime, compressing a gelcauses local and temporal charge inhomogeneity of salt ions. The ions get ”dis-placed” from their equilibrium positions migrate to nearby equilibrium positions.Differential mobility of positive and negative ions produces an electrical current inthe system [61].When the salt concentration becomes significantly lower C ≈ 0.04 M, i.e., De-bye lengths λD ≈ 1.5nm electrostatic interactions become less screened and thesolvent motion becomes increasingly important. In this regime convection of sol-vent molecules carries the ions and produces an electrical potential, that is oftenreferred to as a streaming potential. This has been extensively studied by Prud-nikova and Utz [54].Further lowering the salt concentration (or not adding it at all) and disregard-ing any transient phenomena we arrive at the equilibrium regime that is describedby Donnan type potentials (which will be the primary focus of this work). Thisregime is only observed in the PE gels comprising of a dissociated polyelectrolytenetwork that traps the floating counterions (bearing a charge that is opposite to thebackbone ions). In this regime the potentials arise from the local breakage of theelectroneutrality through the counterion diffusion outside of the gel network. Herethe counterion backbone concentrations and their interactions are the driving forcebehind the build-up of the Donnan type potentials. The counterions appearance isdictated by the pKa/pH of the system, which in turn depends on the local osmoticpressure. The local osmotic pressure could be changed due to the external stimuli4Debye length λD =√εrε0RT2×103F2C0 where C0[M] is the molar concentration, F = eNa =96485.33212Cmol−1 is a Faraday constant, R = 8.31JK/mol is the gas constant. Alternatively,for T = 300K a simple mnemonic formula could be used λD ∼ 0.3[nm]/√c[M]6∆l/l0time, min-∆pH0.00.05 100.40.04Figure 1.2: A schematic representation of the explanation of the piezo-ioniceffect proposed by Sawahata et al. [50, 64] (see text).like applied touch giving rise to the stress-induced ionization of the gel.Stress-induced polyelectrolyte hydrogel ionization and build-up of Donnantype potential was first studied by Sawahata et al. [64], [50]. Sawahata et al. [64],[50] found that an electrical potential could also be generated by a non-uniformcompression of a gel. In his experiment, Sawahata and authors subjected the gelmade of PAA comprising of weak (carboxylic) acid to a uni-axial compression.Upon compression, they observed a reversable change in pH of the gel (Figure 1.2),caused by the appearance of additional protons in the system. As the only source ofappearing protons (H+) was the dissociation of the carboxylic acid, he concludedthat the gel underwent a stress-induced ionization. After subjecting a part of a gelto compression leaving the other part untouched Sawahata et al. [64] observed ageneration of an electric potential. The generation of the potential followed thecompression pattern (as shown in Figure 1.3) linking the stress-induced ionizationwith the generation of the potential difference.7∆l [mm]time, minV [mV]0.0 0.02 410.01.0 VFigure 1.3: A schematic representation of the explanation of the piezo-ioniceffect proposed by Sawahata et al. [50, 64] (see text).The explanation of the stress-induced ionization could be obtained from a purethermodynamic reasoning. In a normal state, backbone molecules remain neu-tral and do not dissociate into positive and negative molecules. However, uponcompression the monomers along the chains become more localized, and, as a re-sult, their entropy decreases. When the entropic incentive becomes strong enough,molecules start dissociating, and the charge density increases. Non-uniform com-pression leads to a non-uniform charge-density, and, as a result, difference in elec-trical potential is generated between two parts of the gel, as shown in Figure 1.4(B),1.3.A missing part in the results of Sawahata et al. [64] was a connection betweenthe ioinzation of the backbone and the gels intrinsic electrostatic potential. Severalgroups [7, 23, 25, 54, 69] dived deeper into this issue by studying a gel immersedin a solvent. They observed an appearance of an electric potential between partsof the gel and the solvent itself. For example, these studies found that gels with anegatively charged backbone (anionic gels) have an intrinsically negative electricalpotential when compared to the surrounding solvent, as shown in Figure 1.4(A).This provided further evidence that the observed potential is an interplay betweenthe diffusion of counterions and build-up of an electrical potential. Unlike Sawa-hata et al. [64] authors connected the observed potential with the Nernst-Donnanpotentials, which is a function of charge density along the polyelectrolyte backbone[54].8Perhaps, the most comprehensive experimental study of the issue was doneby the Utz group, where authors carefully controlled for the ionization and pH ofthe system, and studyied the build-up of an electric potential as well as the os-motic pressure. Prudnikova and Utz [53, 54] performed a series of experimentswhere they formed a gel by co-polymerisation of acrylimide (neutral monomer)and acrylic acid (weak acid capable of dissociation) Poly-(AA-co-AM) in a buffersolution (i.e., keeping the pH constant). By changing the concentration of the AArelative to AM in Poly-(AA-co-AM) they controlled for the system ionization. Byanalyzing the swelling of the gels in water, the authors connected the ionizationof the gel to an increased osmotic pressure. They also observed a difference inelectrical potential and the surrounding solvent, which increased with increasingconcentration of AA (higher ionization). Finally, they studied a compression ofa hydrogel at different ionizations and connected the appearance of the electro-static potential and the applied pressure. These studies further connected the dotsbetween osmotic pressure, applied pressure, ionization and the Nernst-Donnan po-tential [53, 54]The works by [53, 54, 64] provided us with great insights into the behaviourof PE gels in the equilibrium regime. As a result, we observed a connection be-tween the applied pressure, ionization, local osmotic pressure and appearance ofthe Nernst-Donnan potential. However, the nanoscopic understanding behind thelink between the osmotic pressure and ionization of the gel is still missing, as wellas the behaviour of Nernst-Donnan potentials at different strengths of electrostaticinteractions (defined as a Bjerrum length discussed in Section 2.1.2). Experimentsalone are not capable of fully providing this insight. For example, it is a dauntingtask to independently manipulate solvent (and hence the dielectric constant), thetemperature of the solution, and polyelectrolyte ionization in a single experimentalsetup. Coarse-grained molecular dynamics (MD) simulations can be used to fillin the gap. MD simulations represent a bottom-up approach, where the behaviorof each particle in the system is tracked, and collective properties are computed.Most importantly, MD simulations provide a simple framework to manipulate theparameters of the studied system and observe the results.To date, most coarse-grained MD simulations of hydrogels were performed ina setup with an ideal network of polyelectrolyte chains and explicit mobile counte-9Figure 1.4: A: A molecular picture of the Nernst-Donnan equilibrium be-tween two gels with different ionizations studied in experiments. Gelswith a higher degree of negative backbone ionization have an intrinsi-cally lower negative electrical potential when compared to the gels withlower degree of negative backbone ionization. The system consists ofnegative backbone ions (white), floating counterions (green), dissolvedsalt ions (orange and purple) and water molecules. B: A schematicrepresentation of the explanation of the piezo-ionic effect proposed bySawahata et al. [50, 64] (see text).rions floating in an implicit solvent. Usually, the network consisted of a diamond orcubic lattice of cross-linking nodes connected by polyelectrolyte chains [1, 13, 15–18, 33, 43, 55–57, 60, 65–67]. Despite its simplicity, such a setup already revealsinteresting behavior. For instance, Yan and de Pablo [78] demonstrated that in thelimit of weak electrostatic interactions, the gel behaves as an ideal gas as increasingdensity (i.e., compression) of a gel increases the pressure it exhibits. In the limitof strong electrostatics, however, the authors found a somewhat counter-intuitive10Figure 1.5: A: Schematic representation of the interplay between the osmoticpressureΠosm and elastic energy of the networkΠel in the Donnan equi-librium. B: Schematic illustration of counterion condensation. Whenthe Bjerrum length lB of the system becomes bigger than the distancebetween two neighboring backbone macroions a ≤ lB, the counterionsbecome bound to the backbone ions and stop contributing to the os-motic pressure. The system consists of negative backbone ions (white),floating counterions (green).behavior, where the osmotic pressure exhibited by the system decreased upon com-pression. Yan and de Pablo [78] attributed the decrease to counterion condensation.Erbas and Olvera de la Cruz [15, 16] considered a gas of counterions sandwichedbetween two parallel gel slabs. They showed that the electrostatic energy changesupon compression of the gel at fixed strength of electrostatics, which could againonly be explained by counterion condensation. Building upon the results of Yanand de Pablo [78], Mann and coworkers [43] performed an extensive study of gelsof different sizes, ionizations, electrostatic strength, solvent qualities and other pa-rameters. Their comparison between MD results with polyelectrolyte theory iso-lated different contributions to the overall pressure and validated the Manning the-ory of counterion condensation in the limit of straight polymer chains.The works of Yan and de Pablo [78] and Mann [43] provide a comprehensivepicture of a uniform gel under various conditions, while the work of Erbas and11Olvera de la Cruz [15, 16] provides insight into the behavior at the gel-counterioninterface. The question arises how one could combine those studies and create amolecular dynamics setup that provides insight into the piezo-ionic effect.1.2.2 ObjectivesAbove we looked through the body of the available literature and concluded thatthere is a clear gap in understanding of the piezo-ionic effect. In the regime of noadded salt (i.e., no transient effects of migrating salt ions) the relationship betweenionization, Donnan potential as well as the pressure gradient is still unknown. Alsounknown is the dependence of the Donnan potentials on the strength of electrostaticinteractions (i.e., Bjerrum lengths). Based on that we have outlined three mainobjectives for this part of our investigation:• Identify a relationship between the concentration difference (ionization) andthe Donnan potential that is built up at the interface• Identify the change in the Donnan potential at different Bjerrum lengths• Identify the nature of connection of the pressure and potential differencebetween two parts of a sensor1.2.3 Outline of the work on the piezo-ionic effectInspired by the experimental work of Sawahata et al. [64], we consider here agel slab consisting of two parts with different ionization. Unlike the experimentalsetup, where ionization was induced by the lateral pressure applied to one part ofthe gel, the ionization of both parts of the gel in our work is different by design.When this system is allowed to equilibrate, an electric potential difference builds upat the interface, and the corresponding lateral pressure difference can be obtainedor calculated. In Chapter 4 we analyze the ion density, potential and pressureprofiles and show how they vary with temperature and ionization degree. We thenshow that the electrostatic potential for different gels can be collapsed onto a singlemaster curve given by the classic expression for Donnan equilibrium. We finallyshow that the electrostatic potential is proportional to the lateral pressure differencebetween the two gels.12Figure 1.6: Schematic representation of the polyelectrolyte gel diode in Arectifying and B blocking mode. Due to applied voltage, water presentin the gel diode undergoes an electrolysis process generating H+ andOH– ions at the cathode and anode, respectively and, simultaneously,produces electrons in the external circuit. The asymmetry of the poten-tial drop either attracts the counterions to the interface (thus allowingOH– , H+ groups to recombine and rectify the current) or pushes counte-rions to electrodes (thus blocking OH– , H+ groups from recombinationand blocking the current).1.3 Polyelectrolyte Gel Diodes (PGDs)1.3.1 Literature reviewUnlike regular p-n diodes, PGDs remain poorly understood despite being theorizedmore than 70 years ago by Lovrecek et al. [41]. Semi-conductor diodes boast highrectifying ratios, high operating voltages and currents. Polyelectrolyte gel diodes,on the other hand, are still unable to operate in such regimes despite having superior13mechanical properties [26]. Largely, this is due to the lack of understanding and thecomplexity of the electrostatic phenomena occurring inside of the polyelectrolyte.Below we summarize the available body of literature on PGDs as well as outlinethe gaps still in the knowledge yet to be investigated. 5Most studies of PGD performance so far could be divided into three main cate-gories: experimental (mainly I-V current-voltage curve measurements), theoretical(continuum based electrostatic calculations employing Poisson-Boltzmann equa-tions) and molecular simulations (studying the behaviour of individual atoms at amolecular level).The authors of several works [26, 71, 76, 82] have synthesized PGDs with in-triguing chemistry , and have experimentally observed a close connection betweenthe ion structure, electrostatic potential at the electrodes and the performance ofthe PGD. For example, they studied PGDs with variable salt concentration andfound that diodes rectify currents better with less salt present. This implies thatthe performance of the PGD is closely connected with the ions that are native tothe gel (counterions and backbone ions). Authors also observed, that upon applica-tion of constant voltage a steady state current was observed, implying a continuouselectrochemical reaction to happen in the gel powered by the applied voltage. De-spite providing some qualitative evidence on the behaviour of PGD, no quantitativeanalyses and predictions have been made. These gaps were partly addressed viatheoretical models.Despite a growing number of experimental observations, there have been fewtheoretical studies and quantitative calculations of PGDs [42, 70, 77]. One ofthe most comprehensive models was one developed by Yamamoto and Doi [77].Yamamoto and Doi [77] used a simple model based on the water electrolysis,Poisson-Boltzmann and Nernst-Planck equations to analyze the electrostatic po-tentials inside the PGDs to predict the I-V characteristics of PGD. They argue,that in the static regime where there are no electrochemical reactions (hypotheti-cally) at the electrodes the floating counter ions in the gel redistribute themselves5We note that there have been three main types of flexible devices reported in the literature capa-ble of rectifying electric currents. Hayward group [31] synthesized a junction between ionoelastomernetworks with different mobilities of its constituents. Another group used an assymetric oxidationof the electrodes [71]. The most widely used type so far has been one that was based on waterelectrolysis reactions happening at the electrodes, that is the focus of the current investigation.14to minimize the Helmholtz free energy. Thus upon applying a constant potentialat the electrodes the ions migrate from the interface between two gels creating adepletion zone and thus generating large potential drops at the interface. Whenwe allow (hypothetically) the reactions to take place due to water electrolysis H+and OH– are created near the electrode surfaces. These ions then diffuse throughthe gel to the interface between two gels where they recombine and form waterH++OH– −−→H2O and create a current with density J. Due to the inherent asym-metry of the gels, those ions are either heavily screened via floating counterions inthe case of negative bias. In the case of forward bias, however, the screening isdone via the fixed backbone ions (which are much less mobile). Hence, ions freelydiffuse to the interface ”uninterrupted” resulting in an exponential I-V curve. Thescreening (in case of the negative bias) via the floating counterions supresses thecurrent so it scales as J ∼ −√−V , essentially rectifying the current. When thecurrent in both regimes is written-out as a function of charge ionization Yamamotoand Doi [77] make one of the most counter-intuitive claims of the work: with lowerbackbone ionization PGDs will have a higher current in the forward bias, as wellas lower currents in the backward bias regime, i.e., less ionization means betterrectification.The Yamamoto and Doi theory provides a comprehensive model for the work-ing mechanism of a PGD. It is quite remarkable given the relative simplicity ofits assumptions. These assumptions however need to be tested to prove they canpredict the experiment. Namely, there are three blocks of these assumptions: in-puts of the model, the model itself and the approximations used there, as well asthe outputs of the model. For example the Yamamoto model assumes a constantconcentration of the gel backbone ions, which is not true in the case of flexiblenetworks, that can move or collapse. The model itself needs to be also tested. Thisis important because Moy et al. [48] have shown that generally Poisson-Boltzmann(PB) models similar to the one used by Yamamoto and Doi [77] are not univer-sally applicable. These assumptions will be stress-tested via molecular dynamicssimulations. Additionally, potential drops mentioned in the present work couldqualitatively change depending on the strength of the solvent used in the synthesisof the PGD. Hence, a more detailed analysis like molecular simulations is neededto demonstrate the applicability of the PB model as well as to measure the potential15drops at the interface of the gel.There have been limited attempts to tackle PGDs through molecular level sim-ulations. Some studies focused on tangential topics like ionic conductivity of gels[38]. These studies tested PE gels in different operating regimes, e.g., weak andstrong electrostatic coupling (with observed counterion condensation). Addition-ally, Li et al. [38] found that counterion condensation may impede the ion transportin the presence of the electrostatic field, hence implying that classical PB models(like the one used in [77]) will no longer be valid. One work has investigatedPGDs using grand canonical Monte-Carlo simulations [36]. The authors observedthat as a positive voltage is applied (forward bias), small ions start moving acrossthe gel to induce electric current flow. Additionally, with negative voltage (reversebias) applied, the current flow became significantly weaker. Despite being in goodagreement with experimental observations, the methodology employed by the au-thors raises questions about its validity. For example, the authors employed timedependent variables for calculating the electric current, which in the context of aMonte-Carlo simulation has no physical meaning. There remains thus a need for adeeper understanding of PGDs at molecular level.Despite, the nanoscopic insight provided by the MD simulations, some aspectsof the work of PGD are not obtainable in the simulations. Hence, to further stresstest the assumptions made by Yamamoto and Doi [77] more experimental workneeds to be done. For example, the electrochemical reaction times, assumed to beinfinitely fast in Yamamoto and Doi [77], are not. This becomes especially impor-tant when time-dependent voltages are applied to the electrodes of the PGD. An-other aspect of Yamamoto theory, such as their main prediction (I-V) curve haven’tbeen tested systematically, where the predictions of the theory were directly com-pared to the experimental findings. These gaps will be addressed via systematicexperimental study of the PGD controlling for the gel thickness, gel ionization andetc.1.3.2 ObjectivesBased on the literature review we identify that assumptions made in the theoreticalmodelling of the PGDs are yet to be tested. For example, the assumption of a16constant polymer backbone ion concentration throughout the gel is not applicable,since the network even though less mobile than the floating ions is still able tomove. Hence we have the following objectives to address those gaps:Molecular Dynamics investigation• Understand how adequate the PB model is to predict the ion density andelectric field profiles of the PGD– How do density profiles compare between the MD simulation and thePB model?– How valid are the electric field profiles obtained in the PB model?– How do the profiles change for gels with different ionizations and atdifferent Bjerrum lengths• Understand how to modify the PB model to model the PGD more accurately?Experimental study• Understand the basic assumptions of the Yamamoto theory– How is the current created in the system?– How does the system equilibration to the steady state coincide with thediffusion time of the H+ and OH– ions predicted by the Yamamototheory• Stress-test the predictive power of Yamamoto theory by comparing the I-Vcurves to the experimental results1.3.3 Outline of PGD workHere, we expand on coarse grained molecular dynamics (MD) simulations of poly-electrolyte gels [38, 75] to better understand the behaviour of PGDs at nanoscale.The present work combines continuum electrostatic modeling with MD simula-tions for a ”side-by-side” investigation of PGDs in the static regime, i.e. with noelectrolysis happening at the boundaries. By comparing the PB results with thecorresponding MD results, the present work carefully tests the validity of the PB17model used by Yamamoto and Doi [77]. Furthermore, our work develops an im-proved PB model using information obtained via MD simulations. This PB modelcan describe the experimental results more accurately and act as a compass forformulating materials to make better soft electronic devices.Modelling of the PGD can be divided into two steps. The first step is to modela system with no H+ and OH– ions created at the electrodes. At this step, only astatic picture of ion distributions is considered. In a second step, electrolysis andion flow is added into consideration. The ion creation and reconfiguration is mod-elled dynamically, and currents are analyzed as a function of applied voltage. In thepresent work, we focus on the first step and carefully investigate the ion distribu-tion due to electric fields using Poisson-Boltzmann (PB) and Molecular Dynamics(MD) models while not considering electrolysis and currents at the electrodes.By comparing the PB results with the corresponding MD results, the presentwork carefully tests the validity of the PB model used by Yamamoto and Doi [77].Furthermore, our work develops an improved PB model using information obtainedvia MD simulations. This PB model can describe the experimental results moreaccurately and act as a compass for formulating materials to make better soft elec-tronic devices.Experimental investigation, will focus on systematic study of the PGD. Severalsets of PGDs will be created of variable degree of ionization and length. Then thesystem will be subjected to the series of step biases in applied voltage. Observinga current created due to the applied voltage as well as how it scales based on thedegree of ionization and the size of the material will allow us to further-stress testthe Yamamoto theory.1.4 General outline of the thesisThe following part of the thesis outlines fundamental theoretical principles thatgovern the piezo-ionic effect and mechanistic function of PGD. We attempt to un-derstand the Polyelectrolyte gels at molecular level, dive deep into the theory ofpolymer physics, Donnan potentials as well as derive the PB model and introduceimportant parameters like Debye length λD and Bjerrum length lb (Chapter 2).Consequently we describe the methodology behind the current work (Chapter 3).18This will be followed by Chapter 4 where we present the results of the MD simu-lations of soft sensors. In Chapter 5, and 6 the MD and experimental findings onPGD are presented.19Chapter 2Theoretical backgroundIt doesn’t matter how beautiful your theory is, it doesn’t matter howsmart you are. If it doesn’t agree with experiment, it’s wrong.” — Dr.Richard Feynman2.1 Theoretical aspects of polyelectrolyte gels2.1.1 A molecular view of polyelectrolytesA polyelectrolyte hydrogel is a three dimensional network formed by charged poly-mer chains that can swell in water and retain a significant amount of water whilemaintaining their structure [49]. Based on the origins of the chains in the network,the hydrogels could be split into two categories: synthetic or natural.Most widely investigated synthetic PE hydrogels would be, poly (acrylic acid)(PAA), poly (methacrylic acid) (PMAA), and their derivatives. Unlike the natu-ral polymers, the wide variety of monomers enables one to prepare the hydrogelwith desired physical properties for a given application. Furthermore, the sidecarboxylic groups of PAA can be easily modified with various functional groups,such as imide, amine, or attached to other molecules or bioactive agents. This willmake hydrogels sensitive to the external stimuli and make them ”intelligent.” Forexample electrically-responsive PAA polyelectrolyte hydrogels were used as drugdelivery systems [63].20Figure 2.1: An example of a flexible PE. Left: constitution formula for sul-fonated polystyrene and its sodium counterions; right: a coarse-grainedrepresentation that forsakes the chemical detail to focus only on thephysics of the system (taken from [45])Natural hydrogels are ubiquitous in our everyday life. These polymers includepoly-saccharides, such as alginate, heparin, heparan sulfate, hyaluronic acid (HA)and chitosan, proteins/ polypeptides, such as collagen, gelatin and DNA. Both nat-ural and synthetic hydrogels could be further divided into anionic and cationicbased on their pendant acidic groups (e.g., carboxylic acid, sulfonic acid) or basicgroups (e.g., amine groups) on the polymer backbone.Cross-linking is a process that allows a solution of unconnected chains to coa-lesce and form a network. It could be characterized as either chemical or physicalin nature. Chemical or covalent bonds are usually created from a solution of lin-ear polymers. The irradiation of polymer with high energy radiation like UV lightin aqueous solution results in formation of free radicals on the polymer chains.Free radicals initiate the formation of covalent bonds or grafting polymerization,and form a network structure. PAA and other polymers have been cross-linked byradiation method [30]. Alternatively, in some cases, non-covalent cross-linkingreversibly transforms a solution of polymers into a gel. The cross-links of the net-work have a physical origin (hydrogen bonding, Van der Waals forces) and there-fore are sensitive to variations of temperature, pH, ionic content, etc. [39, 40].21Figure 2.2: An example of a rod-like PE. Left: constitution formula forpoly(paraphenylene) with iodine counterions; right: a coarse-grainedrepresentation (taken from [45])2.1.2 Introduction to Polymer PhysicsThe concept of coarse-grainingPolymers are chains of monomers bonded together into long strands via cova-lent bonds. The elementary units of chains called monomers could be any de-sired chemical group. For example binding together CH2 groups can produce verywell known plastic called polyethylene. The chemistry of monomers, turns outto be secondary in determining the chain conformation. Thus, researchers usuallystudy simplified systems of uniform beads tied together into long strands. Havingstripped the system of its complexity, this way, and leaving only the essential ele-ments that will represent the system we will explore below the basics behind thepolymer physics of neutral and charged polymers (polyelectrolytes) that are funda-mental to study the PE hydrogels. The introduction will be minimal and interestedreaders may find more details in introductory books on the subject by Rubinsteinand Colby [58], Flory [21], and De Gennes [11].Physics of Neutral polymersPolymer chains and end-to-end distanceWe will start by studying the simplest version of polymer chains, that are n22Figure 2.3: A schematic representation of an ideal chain (taken from [59])non-interacting spheres interconnected by bond lengths li ≡ |−→r i|. The end-to-enddistance of such a chain is:r2E = (−→rE )2 = (Σni−→li )2 = nl2+2Σ j>i−→li−→l j (2.1)Averaging throughout the ensemble of different conformations, and accountingfor the independence of individual units of the chain (and hence expected value of(−→li ,−→l j ) = l2 < cosθ >= 0) we will find the value for the mean-square end-to-enddistance R2E =< r2E >:R2E =< r2E >= nl2+ l2 < cosθ >= nl2 (2.2)A more realistic model of a polymer is a Freely Jointed Chain (FJC) model,where the chain is modelled as segments (called Kuhn segments) of length b joinedtogether via a joint. The maximum length of such a chain is Rmax =Nmb, where Nmis the number of joined segments. The mean-square end-to-end distance of FJC is:< R2E >= Nmb2 = bRmax. We note that in present work, where we study a highlyidealized and generalized polyelectrolyte model we assume our basic units to befully independent, hence the beads of our polyelectrolytes to be Kuhn segments,i.e., l ∼ b.Tension blobs in ideal chains:To further study polymer chains we will need to introduce a concept of a ten-23Figure 2.4: Left: a schematic representation of thermal blobs in a polymerunder tension; middle: chain extensions for different models of polymerchains; right: comparison of a chain extension between WLC model andexperimental results for λ -DNA (taken from [59])sion blob. A thermal blob is a block of g monomers of a polymer chain whosemotion is fully dictated via a random motion and is not affected by external factors(like stretching of the whole chain). Thus a polymer chain could be seen as a setof thermal blobs (of size ξ ) tightly connected with each other (Figure 2.4).Stretching of a chain won’t perturb the motion of its monomers inside the blob(provided it is small enough). Instead, it will reorganize the blobs themselves intoa straight line. The monomers inside of the blob will continue following an idealchain confirmation (ξ = b2g), as they will remain unperturbed. The number ofblobs could be easily calculated by dividing the number of monomers by the num-ber of monomers per blob Nm/g, and hence the chain length as Rx = ξNm/g.Note, that the loss of entropy due to alignment creates a positive contributionto the Helmholtz free energyF , and acts as a source of tension that the chain willstart exhibiting upon stretching. The free energy of the blob could be calculatedusing the equipartition theorem of kBT per each degree of freedom (i.e., per eachblob) as shown in Equation 2.3.F = kBTNmg= kBT−→RE 2Nmb2(2.3)Real chainsWhile the previous equations considered ideal chains that are modelled asnon-interactive with an external medium the real systems are far from such. Themonomers in the chains interact with themselves as well as with the solvent they24are submerged into. A parameter that characterizes the strength of the interactionsbetween the monomers in the chain and the external medium is the excluded vol-ume vex. Excluded volume is a mathematical expectation of a monomer size basedon the probability of finding other solvent molecules around it. It is used as a mea-sure to quantify the conformations of real chains depending on their interactionsbetween its monomers. Comparing vex with b3 one could observe several scenarioshappening in real systems:vex = b3 In this limit we have an ”athermal” regime where our beads act as hardspheres strongly repulsing each other, independent of temperature (e.g. polystyrenein ethyl benzene)0 < vex < b3 We have a good solvent, where monomers don’t distinguish the attractionwith other monomers or solvent (for example polystyrene immersed in ben-zene)vex = 0 If we further reduce the excluded volume we reach the θ temperature wherechains have nearly ideal conformations−b3 < vex < 0 Further reduction of the excluded volume means that attractive wells domi-nate the interactions and monomers stick close together due to net attraction(e.g. polystyrene in ethanol)vex =−b3 The limiting case of poor solvents is the case of non-solvent (e.g. polystyrenein water)Good solvents:Interactions with a solvent produces an additional interaction energy-borneterm in the Helmholtz free energyFint . Based on that, one could calculate the equi-librium conformations of the chains using the Flory approach. For a good solvent,for every interacting monomer an additional kBT will be added to the Helmholtzfree energy. It represents interaction energy of excluding solvent molecules. Thefraction of the interacting monomers is vexNm/R3E . Given Nm monomers the totalfree energyFtot could be summarized below in:Ftot =Fint +Fent = kBT (vexN2mR3E+R2ENmb2) (2.4)25Figure 2.5: A schematic representation of an ideal chain (top) vs real chain(bottom) of the same counter length Rmax = Nmb stretched by the sameforce−→fx (taken from [59])After minimizing the equation for optimum size RE one finds the Flory radiusRF :RF = v1/4ex b2/5N3/5m ∼ N3/5m ∼ Nνm (2.5)A more generalized form of RE ∼ Nνm has been confirmed via numerous com-puter experiments and has been widely used. It is ν = 1/2 for an ideal solvent, andν = 0.588 for a good solvent [9, 11, 19, 21, 37, 58, 83].Physics of charged polymersFor now we have worked only with neutral polymers with no charge on theirbackbone. General concepts outlined above could also be applied for the weaklycharged chains with: f kBT lbb/ f  kBT . Here we introduced an important quan-tity called the Bjerrum length lb = e24piε0εskBT which is the length scale at which thestrength of the electrostatic interactions is equal to kBT . It reflects the strength ofthe electrostatic interactions relative to the entropic forces present in the system.When chains are weakly ionized the entropic interactions dominate the motion of26the counterions and they could be described as a cloud of an ideal gas-like atoms.As the degree of ionization is increased, electrostatic interactions become increas-ingly important and the ideal-gas approximation is no longer valid.We will now consider the chain stretching mechanism that was considered forideal chains and apply it for the chains with charges in good solvents. Our goalhere is to derive similar equations that will produce the scalings of the chains anddescribe the equilibrium configurations of chains in good solvents. To do that wewill first recall the classic thermal blob picture and subject this chain to a uni-axialtension. We will then introduce the electrostatic interactions in the picture via theconcept of an electrostatic or polyelectrolyte blob. The electrostatic blob is analo-gous to a thermal blob, but applied to a charged system. The competition betweentwo blobs will allow us to define the equilibrium configuration of a charged chainunder tension in good solvent condition.Thermal blobs:Lets consider a chain under a uniform tension with force fE =−∇F and divideit into a group of gT monomers (a thermal blob) of a size ξT where fEξT = kBT .This results in a chain of Nm/gT blobs. By definition this chain of size RE willproduce an elastic pressure:ΠE =− fERE/R3E =− fEξTNmgT=−kBTR3ENmgT(2.6)At the same time, in the limit of weak electrostatic interactions, counterions actas a cloud of non-interacting particles exhibiting only ideal gas like pressure ΠC:ΠC = kBT NCI/V = kBTf NmR3E(2.7)Adding the two pressures together 0 = ΠC +ΠE will produce an equilibriumblob size: gT = 1/ f [6, 43–45]. This in return will produce a scaling relation fornetwork of charged strands RE = f 1−νNmb, in contrast to a neutral chain whereRE = NνmbPolyelectrolyte blobs:As we increase the Bjerrum length, or increase the degree of ionization, Coulom-bic forces become increasingly important and start to dominate the interactions. To27characterize this transition, a new length scale is introduced: polyelectrolyte (PE)blob with gPE monomers and size ξPE = bgνPE . The electrostatic blob remainsinsignificant as long as gPE < gT , which leads to Equation 2.8. The criterion inEquation 2.8 presents a transition between two regimes. The first regime is domi-nated by the entropic forces where the thermal blob is larger than the electrostaticblob gPE < gT ; in the second regime the PE blob becomes larger than the thermalblob and electrostatic interactions become more dominant in the system leading tocounterion condensation discussed below.f > (lb/b)−1/ν (2.8)2.1.3 Counterion condensationAt low temperatures or in the presence of less polar solvents like n-propyl alcohol[2],electrostatic interactions dominate over thermal energies and the classic picture ofion behaviour inside of the gel is no longer valid . This transition is usually quanti-fied by introducing two length scales: the Bjerrum length lb; as well as the distancea between two consecutive charged monomers separated by 1/ f − 1 atoms in thechain [44]. For the case of a good solvent, a could be estimated by a = b f 1−ν 1,where ν = 0.588 and b = 0.97 σ is the average bond length.When a lb, which occurs at high temperatures or in the presence of highlypolar solvents like water, electrostatic forces are weaker than the entropic forces.Following Erbas et al.[14] we will refer to this situation as weak coupling. Bydecreasing a (e.g. increasing the ionization degree) or by reducing lb (e.g. bychoosing a less polar solvent) we can create a system where lb becomes larger thana. In this regime all counterions present in the system become effectively bound totheir corresponding macroions on the backbone in a phenomenon called counterioncondensation[43]. The bound counterions neutralize their respective backbone ionseffectively reducing the fraction of free counterions f → fe f f . Manning’s theoryprovides a simple expression for the effective ionization of the gel as a function ofthe Bjerrum length:1As the length scales are smaller than the electrostatic blob [58]28fe f f = f , for lb ≤ af/(lb/(b f−ν)) , for lb > a, (2.9)It is worth noting that the Bjerrum length lb is hard to vary experimentally asit usually requires changing the dielectric constant via changing the solvent (ex-plicit temperature dependence is negligible as it is usually negated by the implicitdependence of dielectric constant on temperature). Molecular simulations providean easier way to study different scenarios where both parameters a and lb could bevaried. For this reason they have been widely used to perform comprehensive stud-ies of PE systems that would be cumbersome to perform in experimental settings[14, 43, 75].2.1.4 Donnan equilibriumAs was already mentioned, the PE gels comprise of an elastic network of bothcharged macroions (that depends on the degree of ionization f ) and neutral monomers.Every charged macroion is compensated by an itinerant counterion to preserveelectroneutrality. It seems like counterions would favor to leave the hydrogel andescape into the external solvent as this would increase their translational entropy.However, this would break local electroneutrality. Thus, in case of an immovablegel network the counterions remain within the cell but exert an osmotic pressureonto the surroundings which could be calculated via an equation for an ideal gasEquation 2.7, Equation 2.10.Πosmideal = kBTNciV= kBTf NmV, (2.10)If the gel network is stretchable, the osmotic pressure swells the network. Thestretching of the chains decreases their entropy, which creates an elastic pressure asillustrated in Figure 1.5(A). The elastic pressure produces mechanical work ∆Wmechdue to the expansion of the network.Besides the mechanical work, entropic forces could potentially break the electro-neutrality of the system but only on small length scales l < Ree, where Ree is theaverage distance between any two cross-linking nodes. The breaking of electro-29neutrality produces an electrical potential ∆ϕ , also called as Nernst-Donnan po-tential, which follows readily from equating the electrochemical potentials of thecounterions of charge q inside and outside of the gel µin− µout = 0, where µ =µo+ kBT ln(c)+qϕ+µexcess:kBT ln(co/ci) = q∆ϕ+∆µexcess, (2.11)where ci and co are the average counterion concentrations inside and outside ofthe gel and µexcess is the excess chemical potential due to excluded volume interac-tions present in the system due to solvation.2.2 Theory of Polyelectrolyte Gel Diodes2.2.1 Poisson Boltzmann equationTo model the PGD consisting of a positive and negative gel sandwiched betweenthe electrodes, it is sufficient to consider only half of the system due to the inherentsymmetry. On the right side of the interface we have a positively charged backbonewith a concentration of the fixed charges c f and itinerant positive and negativecounterions with a concentration of c0. Assuming a Boltzmann distribution of theitinerant ions and a uniform distribution of the fixed ions, the electrostatic potentialφ(z) 2 along the (z-)axis of the diode can be described asd2ϕ(z)dz2=−(ρit(z)+ρbb(z)) · ( qεdielkBT ) == 8pilbc0 sinh(ϕ(z))−4pilbc f .(2.12)The concentration c0 is defined so that∫ z=bz=a dzc0 exp(ϕ(z)) = c f L/2, where L isthe length of the simulation domain and ϕ is measured in dimensionless units (i.e.,ϕ = qVkBT ). By introducing the inverse Gouy-Chapman length l−2GC = k20 = 8pilbc0, aswell as k2f = 4pilbc f , where lb = kCoulq2/εdielkBT is the Bjerrum length with εdielthe dielectric permittivity of the medium and kCoul the Coulombic constant, the PBequation can be rewritten as2This is a dimensionless potential in units of kBT .30d2ϕ(z)dz2= k2o sinh(ϕ(z))− k2f (2.13)Note that the parameters c0,k0 are not known a priori unlike c f ,k f , hence tosolve the equation above an iterative approach must be employed. This equationwill be referred to as a standard Poisson-Boltzmann model 3.2.2.2 Transport equationsTwo integral parts of a working PGD are the asymmetry of its structure, as well asappearance and subsequent transport of H+ and OH– ions from the electrodes tothe interface. The asymmetry of the PGD brings an asymmetric electric response ofthe ions inside of the gel, following the PB equation. The transport of H+ and OH–ions, however, is described via Nernst-Planck like equations, that in the equilibrium(i.e., steady state) could be described in Equation 2.14. The Equation 2.14 is acombination of diffusion transport due to inhomogenuity in ion concentration aswell as the drift due to the applied electrical field E =−∂ϕ(z)/∂ z.JH =−DHhyd[∂∂ z cH(z)+ cH(z)∂ϕ(z)∂ z], (2.14)JOH =−DOHhyd[∂∂ z cOH(z)− cOH(z) ∂ϕ(z)∂ z]Where JH ,JOH are the current densities, DHhyd ,DOHhyd - are the diffusion con-stants, of the H+ and OH– ions respectively. Theoretically, these diffusion con-stants strongly depend on the density of the ionic groups in the PE. However, in thelimit of infinte dilution those constants are well known DHhyd = 1.27 · 10−4cm2/s,DOHhyd = 7.15 ·10−5cm2/s 4.3Note that in SI units the Equation 2.12 could be described asd2ϕ/dx2 = (c0e/εε0) · (eeϕ(x)kBT − e−eϕ(x)kBT )− c f e/εε04In regular SI units the equation could be J =−D[∇c+ zFRT c(∇φ)], where F = eNa is the Faradayconstant312.2.3 Yamamoto and Doi theory of PE diodesThis theory provides a thorough explanation on the operation of PGDs. The mainargument lies in the electrochemical reactions at the electrodes where water is splitinto two ions H+ and OH– and their subsequent transport to the interface betweenthe gels. To explain the operation of the PGD we will first look at the static regimewhen the electrochemical reactions at the electrodes are suppressed, followed bythe regime of reverse bias (blocking regime of the PGD) and, finally, the regime ofthe forward bias.Static regime (with currents suppressed at the electrodes)In the static regime, the Yamamoto model treats the ionic groups of the PE back-bone as a uniform distribution of fixed electric charges. The floating counterionsthat have appeared due to the dissociation of the PE backbone, remain close to theirrespective backbone to preserve the local electroneutrality. They re-equilibratethemselves and follow a smooth Boltzmann like distribution. Their profiles couldbe calculated by solving a Poisson-Boltzmann equation Equation 2.12 given that∫ z=bz=a dzc0 exp(ϕ(z)) = c f L/2. The resulting profiles are shown in Figure 2.6.The Figure 2.6 shows the electrostatic potential and ion density profiles of thehalf of the system. It identifies three distinct parts of the diode namely: the counte-rion depletion area at the junction of oppositely charged gels 0 < z < zd , PGD bulkregion (zd < z < ze) and the surfaces of electrodes ze < z < L/2.In the bulk region, as could be expected, floating counterions redistribute them-selves until they fully screen the interactions of the fixed ions and so that c−(z) =c f . Redistribution continues until the potential in the whole bulk region of the PGDis constant (and equal to ϕ f = asinh(c f2c0)). Due to the inherent anti-symmetry ofthe diode, the ion distribution on the left side is opposite to the right side, withpotential −ϕ f generated in the bulk of the left hand side.When a reverse bias is applied −V0 > 0, the electrode on the right hand sidegets a positive charge. To neutralize the charge on the electrode and minimize theHelmholtz free energy counterions from the interface travel to the electrode andform an electric double layer. Here, despite a large applied voltage −V0 >> 1 thepotential felt at the double layer is relatively small (−V0−ϕ f )∼ 1. However, since32neutralization of the electrode surfaces occur at the expense of the ions present atthe gel-gel interface the number of counterions at the junction decreases. Thisresults in large potential drops that are felt at the interface between the gels andeven smaller overpotential at the electrodes.An opposite situation is observed, when a positive bias is applied −V0 < 0.Here the electrode on the right hand side gets a negative charge. The negativelycharged electrode drives the negative counterions away from the surface of theelectrode to the interface between the gels. This creates a large potential drop”that is felt” at the surface of the electrode, populates the ”depletion zone” (at theinterface between the gels) with counterions and lowers the bulk potential ϕ f . Themagnitude of the overpotentials that are ”felt” at the electrode surfaces becomesfundamentally important when the PGDs are studied in the non-static regime.Reverse bias V0 > 0 (non-static regime)Besides the basic electrostatics occuring in the gel, Yamamoto and Doi [77] arguethere is an appearance of additional charges at the electrodes due to electrolysisof water molecules. When a voltage applied to the electrodes is higher than thestandard electrode potentials of water dissociation (1.23 V) a chemical reaction2H2O −−→ 4H+ +O2 + 4e– takes place. The rate of H+ generation follows theButler-Volmer kinetics with JH = J0(e−∆V −1).As was mentioned above, the overpotentials ”felt” at the electrodes (∆V ) areinsignificant in the reverse bias regime. When additional H+ ions appear at theelectrodes even more counterions leave the depletion zone to neutralize the appear-ance of the H+ ions. This results in even larger potential drops that are felt at theinterface (ϕ f = asinh(c f2c0) increases as c0 decreases) between the gels and evensmaller overpotential at the electrodes. Small overpotentials result in insignificantcurrents JH = J0(e−∆V −1).Moreover, H+ that migrate to the interface between the gels to associate withOH– (we remind here that a similar reactions take place on the opposite electrode).The migration however is heavily screened by mobile floating counterions. Thismeans that H+ don’t change the character of the electrostatic profiles in the system.However, as we have seen above, they do change number of counterions in the de-33pletion layer, and hence the potential drops at the interface which have a profoundeffect at the values of the reverse currents.The migration of H+ is a free diffusion from the electrode to the interface withshort-range repulsion. Its diffusion length defined as λt = −Dhydct/JH > L (verylarge for small currents), where ct =√c2f +4c20. When the system reaches thesteady state, the hydrogen ions become present in the whole system followingthe density profile shown in Equation 2.15. This profile is obtained via solvinga simplified PB model Equation 2.12, combined with the Nernst-Planck transportequations (Equation 2.14) and boundary condition cH(zd) = 0.cH(z)ct=√1+z− zdλt−1 (2.15)To explain the dependence of the currents on the applied voltage Yamamotoand Doi [77] provide a simple scaling reasoning showing that the currents are fol-lowing the relation shown in Equation 2.16. The reasoning is based on calculatingthe number of counterions that are mobilized from the depletion layer to neutral-ize the flock of H+ ions and connecting it to the potential drops at the interface(since that is where it is ”felt” in the regime of reverse bias). Following the Gauss’theorem this number of charges is equal to ec f zd = σ , where σ surface chargehas been applied to the electrodes. On the other hand, ϕ(z) =∫ zd0 E(z)dz, whereE = eσ/(εT ). Hence ϕ f =∫ zd0 dze2c f /((εT )) ≈ 0.5k2f z2d . For larger values ofreverse voltages, overpotentials are insignificant ∆V = ϕ f −V0/2 ≈ 0, hence thepotential drops could be connected to the applied voltage: ϕ f ≈−V0/2.To add currents into the consideration we will now calculate the number ofH+ ions present in the system. For large diffusion lengths λt , local density ofions could be approximated via Equation 2.15: nH(z) ≈ ctz/λt , hence the totalnumber of hydrogen ions in the system is NH ≈ ctL2/(8λt). As was discussedabove, the hydrogen ions travelling from the electrodes to the junction are screenedvia counterions. These counterions arrive from the depletion area that is zd thick,hence NH ≈ c f zd . The depletion of the junction between the gels results in c f >>c0, and hence λ f ≈ λt . This further simplifies the equation for the thickness of thedepletion zone zd = L2/8λ f . Finally, equating two equations for the thickness ofthe depletion zone (i.e.,−V0/2≈ ϕ f ≈ 0.5k2f z2d and zd = L2/8λ f ) gives us a simple34relationship connecting current and applied voltage in the regime of reverse biasEquation 2.16.JHc f Dhyd≈−8√−V0/(k f L2) (2.16)More general form of the potential drops as a function of the number of hydro-gen ions is presented in Equation 2.17.k f L2tanh(ϕ f )≈ k f (L/2− zd)− tanh3/2ϕ f −EL/2k ftanh2ϕ f +k f NHcttanhϕ f (2.17)In the limit of relatively large reverse voltages Equation 2.17 could be approx-imated by Equation 2.18.k f zd ≈−EL/2k f+k f L28λ f−1 (2.18)For the case of relatively small reverse voltage, the bulk electrostatic potentialshave an approximate form shown in Equation 2.19.ϕ f = 1/2log(k f L)−1/2log(√log(k f L)−1+EL/2k f− k f L28λ f(2.19)The solution of the Equation 2.17 is shown in Figure 2.7. Unlike in the for-ward bias regime (discussed below) where currents raise exponentially, reversebias regime produces much smaller currents that follow ∼ √−V0 scaling (Equa-tion 2.17, Equation 2.16).Forward bias V0 > 0 (non-static regime)Having understood the system behaviour when negative bias is applied we nowsubject the system to the positive bias, −V0 < 0. In this regime, protons areappearing (2H2O −−→ 4H+ +O2 + 4e– ) on the left side of the diode, but thenegative hydroxide groups OH– are appearing on the right side of the diode via2H2O−−→ H2+2OH– −2e– with the standard electrode potential of 0.83 V.35Counterion density profiles, as well as the potential profiles remain similar tothe static regime where the electrochemical reactions were suppressed. The OH–still freely diffuses from the electrode to the interface with short-range repulsionwith second virial coefficient c−1t being screened by the opposite charges. Thediffusion length λt =−Dhydct/JOH > L remains large as in the reverse bias regime.However, there are no floating counterions to be able to screen them effectively (asthey bear the same charge), hence this function is being performed by the positivelycharged fixed ions (at least for small voltages). This will have a profound effect onthe I(V) curves as will be discussed below.When the system reaches steady state, the hydroxide ions become presentthroughout the right side of the diode following the density profile described inEquation 2.20. This profile is obtained via solving a simplified PB model, com-bined with the Nernst-Planck transport equations and boundary condition cOH(zd)=zd/λt .cOH(z)ct=√1+2zλt−1 (2.20)Due to large overpotentials (discussed in Section 2.2.3), electric currents inPGDs increase exponentially with an increase of applied voltages for voltages, seeFigure 2.7. For voltages smaller than the bulk potentials V0 < 2ϕ f ≈ 10 [T/e] ≈0.25 V the currents are independent of k fλ0 (i.e., ionization of the gel). This hap-pens because (for small voltages) migration of OH– ions is done via the positivelycharged fixed ions which are highly immobile.In contrast, when the applied bias becomes larger V0 > 2ϕ f ≈ 10 [T/e] ≈0.25 V , electrostatic incentive for the positively charged counterions from the leftside of the diode becomes strong enough to ”pull” them to the right side of diode.After migrating to the right hand side of the diode they screen the hydroxide ions,as well as leave more unmatched fixed ions to screen the migration of H+ (that areproduced on the opposite side of the diode). Now, however when hydroxide ionsbecome screened, the currents are a function of the number of ”migrant” counteri-ons and hence depends on k fλ0. A good approximation of the current dependancein the reverse bias regime is shown in Equation 2.21.36k f L2tanhϕ f ≈ k f (ze− zd)−2tanh3/2ϕ f + k f Nele2c0cosϕ f −k f NOHcttanhϕ f (2.21)For large electric currents the Equation 2.21 could be simplified (Equation 2.22).ϕ f ≈ tanh−1(λ fL[(1+3L2λ f)2/3−1])(2.22)37Figure 2.6: Electrostatic potentials in the equilibrium. Polyelectrolyte geldiodes (PGD) are the double layers of two oppositely charged poly-electrolyte gels that are sandwiched by two symmetric electrodes. Ex-periments have shown that PGDs show very small electric currents forreverse applied voltages (a), V0 < 0, and show very large electric cur-rents for forward applied voltages (b), V0 > 0, for the cases in whichsalt impurity is eliminated [10, 82]. The bottom figures schematicallyshow electrostatic potentials ϕ(z) (the black solid curve) and the lo-cal densities of positive and negative counterions, c+(z) and c−(z), (thered and blue solid curves) in a half of the PGD, z > 0, for the cases inwhich electrochemical reactions at the electrodes are suppressed. Theionic groups of polyelectrolyte gels are treated as uniform fixed electriccharges (these charges are shown by + and - with filled circles) and thedistributions of their counterions are functions of electrostatic potentials(these counterions are shown by + and - with open circles). PGDs showdepletion layers, where counterions are depleted, at the interface regionbetween the two oppositely charged gels, 0 < z < zd , and the surfacesof the electrodes, ze < z < L/2. Electrostatic potentials ϕ(z) are ap-proximately constant cf in the bulk region of polyelectrolyte gels. Bulkelectrostatic potentials cf are (halves of) potential drops across the de-pletion layer at the interface between the two gels and thus are generatedby the fixed charges of the polyelectrolyte gel in this layer (taken fromYamamoto and Doi [77])38Figure 2.7: Voltammogram (I(V)) profiles of PGD in the system without ad-ditional salt. Electric current densities Jele = (= eJH =−eJOH) (that arerescaled by exchange current density eJ0, see Methods) are calculatedas functions of applied voltages V0 (in the unit of the thermal voltage)for reverse (a) and forward voltages (b). The number of hydrogen andhydroxide ions that are accumulated in the bulk of the gel increases withdecreasing the values of k f 0 (that is defined by equation (2)). We usedk f 0 = 1.0 · 106 (black solid curve), 5.0 · 106 (blue) and 1.0 · 107 (red)for the calculations. The value of (rescaled) gel thickness k f L is fixedto 1.0 · 104. The equations that are used for the calculations are shownin Section 2.2.3, 2.2.3. We used a simple scaling argument (shown inMethods) to predict the (black) broken curve in a. The (black) brokencurve in b shows the cases in which all of the applied voltages are dis-tributed to the surfaces of the electrodes (taken from Yamamoto and Doi[77])39Chapter 3Methodology and materials3.1 Method of Molecular Dynamics (MD)3.1.1 Introduction: What is MD?Despite the growing sophistication of experimental techniques some phenomenarelated to motion of atoms on nanoscale remain hard to reach. To fill in the gaps inour understanding Molecular Dynamics simulations were introduced in early 1960sand have become an increasingly powerful instrument in the scientific toolkit.Molecular Dynamics simulations, as Nobel Prize winning physicist RichardFeynman put it, ”explain the motion of the matter through wiggling and jiggling ofmolecules.” They solve the classical equations of motion for molecules modeled asa system of spheres interacting according to a potential energy force field. Differenttypes of spheres describe the different types of atoms found in the molecule. Theresult of a simulation is a trajectory that specifies how the positions and velocitiesof the atoms in the system vary with time.MD simulations provide a link between physical theories described by equa-tions and real systems studied in experimental settings. They are useful in timeevolution studies of soft-matter materials, as they are able to predict configurationsof materials, chemical reactions, conformational changes in molecular behaviour,thermodynamic parameters like Gibbs/Helmholtz free energies.Simulations can examine all states accessed by a molecule at a temporal and40Figure 3.1: The length and time scales of MD simulations compared to avail-able experimental techniques and other theoretical methods (courtesy ofKochmann research group [32]).spatial resolution that cannot be observed experimentally or through a continuumtheory, as shown in Figure 3.1. For example, MD simulations can reach the fastphenomena like formation of electrical double-layer at the boundary between anelectrode and an electrolyte. These phenomena happen in τDL ∼ 100 ns and areoften unreachable for majority of experimental methods1.Molecular Dynamics simulations also allow us to reach the lengths scales thatare cumbersome to reach in experimental settings, and are fundamentally neglectedin the continuum theories. These length scales are on the orders of the Debye lengthλD ∼ 1 nm, and are important for studying the interface phenomena in Polyelec-trolyte Gel Diode, or providing potential profile in the piezo-ionic sensors.Below we will dive deeper into the core of any MD simulations - integratingthe equations of motion, we will discuss how the systems studied in MD simula-tions are created, how interatomic forces are calculated and provide details on the”nuts and bolts” of the numerical simulation. The introduction will be brief, and1According to Maxwell-Wagner, this time scale could be calculated from the Debye length of anelectrolyte and the corresponding diffusion time: τDL = λ 2D/D ∼ 100 ns, where D is the diffusionconstant, and λD is the Debye length.41readers can find more details in one of the most comprehensive books on MolecularDynamics by Allen et al. [4].3.1.2 Core of MD modeling: equations of motionUnlike the regular modelling theories, where typically one has a single elaborateequation to solve and predict the behaviour of the material, MD simulations workwith a large set of simple Newton’s equations of motion (Equation 3.1). Theseequations are applied for all the particles in the system and represent a coupled sys-tem of N (often millions) 2nd order Ordinary Differential Equations. The systemof equations is then integrated over time to compute the trajectory of the system.mid2ridt2= Fi(r1,r2, . . . ,rN) (3.1)Here ri are the position vectors and Fi are the forces acting upon the N atoms in thesystem. The forces are derived from potential energy functions, U(r1,r2, . . . ,rN),representing the potential energy of the system for the specific geometric arrange-ment of the atoms:Fi(r1,r2, . . . ,rN) =−∇riU(r1,r2, . . . ,rN)The algorithm of performing an MD simulation is illustrated in Figure 3.2 andcould be described in three major steps:• Create initial configurations: where we need to setup initial coordinates ofevery atom in the system• Calculate forces: where we calculate the forces applied to every atom, bysetting the interaction ”strength” through inter-atomic potential• Integrate equations of motion (Equation 3.1) for every atom for every timestepin the trajectoryNewton’s equations of motion conserve total energy, linear and angular mo-mentum of the system. In statistical mechanics this system is characterized as onethat follows the microcanonical ensemble (NVE). The ensemble of integration will42Figure 3.2: The algorithm of an MD simulation. First, initial configurationsare created for every atom in the system. Then forces applied to ev-ery atom are calculated by setting the interaction ”strength” throughinter-atomic potential. Then equations of motion (Equation 3.1) are in-tegrated.become very important when integrating constraints and thermostats will be dis-cussed in the next subsection (Section 3.1.3).3.1.3 The ”nuts and bolts” of the MD simulationSimulation workflowTypically, the coordinates of an initial configuration of the system can be eitherobtained from the experimental results, or designed to represent the physics ofthe system. The initial positions of the atoms have to be carefully set to avoid anyoverlaps in the atom positions, that can potentially generate large interacting forces”blowing up” the system causing the simulation to fail.The next step in the workflow is to equilibrate the atoms in the system. Thepurpose of the equilibration phase is to allow the system to evolve from the startingconfiguration to reach equilibrium. This is done by assigning initial velocitiesfrom a Maxwell-Boltzmann distribution at a low temperature to each atom and43integrating Newton’s equations Equation 3.1 to advance the system in time. Newvelocities are assigned at each step until the desired temperature is reached. Inparallel with temperature, other parameters like structure, temperature, energy arefollowed. By analyzing the dispersion of the observed parameters we measure ifthe system is in steady state. The system will have reached the steady state whenthe measurement error of the sampled observables of these properties are within anacceptable range.During the production phase, the simulation is run for the desired time length(ps to ns) and trajectories (atom positions and velocities) collected. The trajectoriescan then be used to compute thermodynamic quantities (temperature, energy) andvisualize conformational changes. The atomic positions could be averaged andstatic potential and density profiles could be generated.Numerical methods behind integrating equations of motionNumerical integration of Equation 3.1 will advance the system by a small timestep δ t. During this time step the force on the system is considered approximatelyconstant. If δ t is small enough, the solution will be a reasonable approximation.There are many algorithms for integrating the equations of motion using finitedifferences: like Euler method, RK5 and etc [4]. The Verlet algorithm is the mostwidely used method in MD for integrating the equations of motion. The Verletalgorithm uses the positions and accelerations at the current step together with thepositions from the previous step, to calculate the new positions for the next stepEquation 3.2 is used.r(t+δ t) = 2r(t)− r(t−δ t)+δ t2a(t), (3.2)v(t+12δ t) =1δ t(r(t+δ t)− r(t))The Verlet algorithm is easy to implement and requires storage of only two setsof positions r(t) and r(t− δ t) and the acceleration a(t). A variation on the Verletalgorithm, is so called leap-frog algorithm (Equation 3.3).44r(t+δ t) = r(t)+δ tv(t+12δ t), (3.3)v(t+12δ t) = v(t− 12δ t)+δ ta(t)One disadvantage of the leap-frog algorithm is that the positions and velocities arenot synchronized and so it is not possible to calculate kinetic energy. Another formof the Verlet algorithm, called velocity Verlet, can compute positions, velocitiesand accelerations at the same time.There are no well-defined rules for choosing the time step in MD simulations.A time step that is too small will be prohibitively long to run. Too large a time stepmay cause instabilities in the integration algorithm due to high energy interactionsproduced from overlapping atoms. As a rule of thumb, the time step should be anorder of magnitude less than the fastest motion in the system. For polymers, thisis the C−H bond stretching at around 10 fs. This gives an integration time step of1 fs which explains why MD can only access nanosecond timescales if atomisticsimulations are performed.Coarse-grainingTo increase the timescales accessible by computer simulations one can forget aboutthe chemical complexity of the material provided they are not of interest in thesimulation. Thus C−H bonds or even the whole C−H2 units could be representedas a single bead. This significantly reduces the computational cost of the system,as well as allows us to run simulations much longer that are not constrained bythe frequency of the C−H oscillations. One such model is described below inFigure 3.3.UnitsHaving stripped the system of its chemical complexity leaving only the essentialelements that will represent the physics of the system it is natural to move awayfrom the real units of fs and nm to generalized units used in coarse-grained sim-ulations. One way to do that is to use the Lennard-Jones units derived from the45Figure 3.3: The coarse-graining methodology where the whole C−H3 unitsare represented as a single bead.Lennard-Jones potential (discussed below in greater detail). Instead the energy andlength units are defined by the Lennard-Jones potential values εLJ and σLJ . TheBoltzmann constant kB and the mass of each monomer are set to one. This fixesthe dimensions for other units and are referred to as Lennard-Jones reduced units.It is possible to map LJ units to any polyelectrolyte as the simulation units arenothing more than a rescaled variant of the SI system. To obtain the correct map-ping of the parameters we will use PAA gel as our reference. We will start withdefining a length scale [L] = 1 σ = 1.2 nm. The bead of 1 σ in LJ units representsa single Kuhn segment of PAA. These beads are connected via a bond of 0.97 σ ,which shouldn’t be confused with the covalent bonds present in real systems. Theenergy units come naturally if use the definition of soft-matter as a material where46the interaction energies are of the orders of magnitude of kBT . So for room tem-perature 300 K, our unit of energy becomes [E] = 1 ε = kB ·300K ≈ 4.14 ·1021 J.Since our simulations are operating only with static quantities (in steady state andequilibrium), the units of mass will only correspond to how fast our system will beequilibrated and define our simulation timestep. Here we could use the mass of car-boxylic groups [M] =MC3O2H +3 = 3 ·12 µ+2 ·16 µ+3 ·1 µ ≈ 71 µ = 1.2 ·10−25kg.To hold, the unit of time has to be defined so that [E] = [M]([L]/[t])2. From that weget the missing unit of time [t] = ([M][L]2/E)1/2 = 7.73 ps. The last unit we need isthe one of charge. We choose it to be the elementary charge [q] = e= 1.60 ·1019 C.From that one can calculate any missing parameters and convert them into the realunits, e.g., force [F ] = [E]/[L], or bond elastic constant kFENE = [F ]/[L] = [E]/[L]2.A table of conversion is shown in Table 3.1.Parameter Simulation Value SI valueThermal energy 1 ε 4.14·1021 JDistance 1 σ 1.2 nmTime 1 τ 7.73 psBond elasticity constant kFENE 10 ε/σ2 28.75 pN/nmBond constant rFENE 1.5 σ 1.8 nmTable 3.1: Table of conversion between real and LJ unitsInteratomic potentialsAs already demonstrated it is necessary to be able to calculate the forces of a givenparticle configuration in order to apply any sort of time integration scheme. Thesystems that are investigated in this work use only pairwise additive bonded andnon-bonded forces. Both of them can be calculated as the gradient to a radiallysymmetric potential F =−∇U(r).The Lennard-Jones (LJ) potential was originally aimed to model the interac-tions of noble gases like Helium, Neon, Argon etc. As such, it shows a smallattractive region at larger distances scaling as ∼ r−6. At distances smaller thanthe diameter of the considered “atoms” it shows a strong repulsion, which to avoidcomputational cost was modelled as ∼ (r−6)2 ∼ r−12, as it originated at a time47Figure 3.4: The potential profile for MD simulation.where computers just emerged and calculation power was scarce.ULJ(ri j) = 4ε[(σri j)12−(σri j)6](3.4)here ε is a measure for the strength of the attraction, i.e. the depth of the attractiveregion and σ defines the spatial extension of the potential.Closely related to the LJ potential is the Weeks-Chandler-Anderson potential(WCA) which will be used in the present work. The WCA potential cuts off thepotential right at the minimum rcut = 21/6σ and shifts the curve such that it reacheszero at that point. This potential is then only repulsive, but it smoothly goes to zeroas shown in Equation 3.5.ULJ(ri j) =4ε[(σri j)12−(σri j)6−U(rc)],for r < rc0 ,for r ≥ rc(3.5)where rc = 21/6σ . The LJ parameters εLJ and σLJ correspond to the depth of the48potential well and the distance at which the unshifted potential is zero, respectively.In our work both parameters (εLJ and σLJ ) are identical for all species present inthe system (backbone monomers and itinerant counterions) and are used as units ofenergy and distance. This modified potential was shown to represent the conditionsof a good solvent and is widely used for coarse-grained modelling [43, 45, 46].The Finite-Extensible-Nonlinear-Elastic (FENE) Potential on the other handaims to model a finitely extensible bond between two distinct particles. The FENEpotential reads:UFENE(r) =−12k f r2FENE ln(1−(rrFENE)2), for r < r f∞ , for r ≥ r f(3.6)with kFENE = 10 ε/σ2, rFENE = 1.5 σ resulting in a mean bond length of < b >=0.97 σ ∼ σ . Here kFENE defines the stiffness of the potential, rFENE is the max-imum extension before the bond breaks and r is a shift from the position of theabsolute minimum of the potential.Ions interact additionally via the Coulombic potentialUCoul(ri j) = kCoul−q2εdielri j=−kBT lb/ri j, (3.7)where kCoul - is a Coulomb constant (1 in LJ units), q - the particle charge, and εdiel- dielectric permitivity of the medium. An electric field could applied by addingan additional force F = qE to each charged atom in the system (shown in Equa-tion 3.1) and solving updated equations of motion. To reduce the computationalcost, solvent molecules are treated implicitly through the dielectric constant εdiel[27]. The coarse-grained model described above has been widely used in moleculardynamics and Monte-Carlo simulations by the MD community [14, 16, 79].Integrating constraints: boundary conditionsWe use an orthogonal simulation cell - its initial dimensions are cubic but canchange via the application of a barostat or applied deformation. In some cases,49Figure 3.5: Schematic representation of the periodic boundary conditions,where the systems boundaries are shrink wrapped to avoid the finitesize effects.periodic boundary conditions are employed and the polymer chains can freely crossthem, modelling a bulk material. A minimum image convention ensures monomersinteract with the closest image of the remaining system, see Figure 3.5.Integrating constraints: MD thermostats and barostatsAs was mentioned above integrating Newton’s classical equations of motion asabove involves the microcanonical ensemble (NVE) meaning that in the run of thetrajectory number of particles, volume and total energy of the particles will bepreserved. In some cases however, a more realistic system needs to be studied,where for example temperature and pressure is kept at certain value like 300 K and1 bar. By altering the equations of motion Equation 3.1 we can adjust the velocityof atoms or the dimensions of the simulation box. In the first case, we will create50the so called ”thermostat”. In the second case, a barostat.Thermostatting:Note, that unlike the experimental thermostats, which ”dunk” the system into a”heat bath” the temperature coupling in MD is local. The heat is controlled by ad-justing the velocity on a particle-by-particle basis, rather than allowing heat diffu-sion from boundaries (as would be the case for an experimental set up). However,if the timescale for heat diffusion is considered much shorter than the simulatedtime scales both a local and a boundary coupling to a heat bath become equiva-lent. One way of thermostatting the system is to use the Nose´-Hoover thermostat.Here the equations of motion (Equation 3.1) are modified with an additional termξ which enters into the equation of motion as the strength of a ‘frictional force’ onthe particle so that the acceleration is:ma =−δUδr−mξv (3.8)The additional force comes from a ”fictional heat bath” that creates a heat flowfrom the bath to the system shown in :dξdt= 1/Q(−3kbT/N) (3.9)where is total kinetic energy N is the total number of atoms and Q can be interpretedas a ‘heat bath mass’ and is of the form Q= 6NkbTτ2T , where τT is the ‘temperaturerelaxation time’. It could be shown that a system subjected to the Nose-Hooverthermostat will follow a so called canonical or NVT ensemble (where number ofparticles, volume and temperature of the system is preserved). Unlike the NVEensemble where only the energy of the system is conserved, here the energy of theheat bath and the system is conserved together.Barostatting:Nose Hoover barostat is analogous to the Nose-Hoover thermostat for the isothermal-isobaric ensemble (NPT), where number of particles, pressure and temperature ofthe system is controlled. Now, instead of adding just one degree of freedom tothe system for temperature, an additional degree is added for pressure. A simu-lation package LAMMPS uses a modified version of this barostat to account for51the box changing in shape Equation 3.10. The particle positions are rescaled toset the pressure and this rescaling can be coupled to particular spatial dimensions(i = x,y,z).dLidt= νiαLi (3.10)where Li is the length of the i dimension of the system, νi is an ‘effective strainrate’. The equation of motion for position is modified to incorporatedridt= vi+νi(ri−Rcm,i), (3.11)ma = −δUδr−mξvwhere Rcm,i, i is the i’th Cartesian component of the center of mass of the system.The variable ν has the following dynamicsdνidt=V (t)t2BNkbT(Πii−P)Here, Πi j are the elements of the total pressure tensor, P is the fixed external pres-sure, and tB is the pressure relaxation time.A note on used softwareMolecular dynamics runs were performed using the LAMMPS code [51]. The MDtrajectories were visualized using the VMD code [29] and analyzed using the MD-Analysis package [24, 47]. The code base employed in the present work is sum-marized in the PyGels library [74]. All results are reported in reduced simulationunits.3.1.4 Details of implementationPiezo-ionic effectOur simulation setup comprises of a box of 18x4x4 unit cells forming a cubiclattice of crosslinking nodes. The unit cells have lattice constant l = 50σ andare connected by polyelectrolyte chains with Nm = 100 monomers. Due to the52translational symmetry of the system there are three polymer chains per each cellwith monomer number density c = 3 · 100/503 = 2.4 · 10−3σ−3. A fraction f ofthe monomers in each chain are ionized, i.e. they bear a negative charge q = −1.Monomers on the left side of the box are ionized with the ionization degree f1,and monomers on the right side with the degree f2 > f1. After the network iscreated, Nci counterions are inserted into the box at random positions to achieveelectroneutrality. Monomers in the network are tethered to static nodes at the leftand right boundaries of the simulation box. Periodic boundary conditions (PBC)are employed in lateral directions only. Reflective walls are put at z = 0 and z = Lzto prevent atoms from escaping the box. Several snapshots illustrating the systemare shown in Figure 4.1 2.To further optimize the computational performance, the electrostatic cutoff dis-tance was set to relc = 17σ , above which longer-range electrostatic interactions arecalculated via a Particle-Particle-Particle Mesh (PPPM) Ewald solver for slab sys-tems [80]. All simulations were performed in the NV T ensemble, with a timestepof integration dt = 0.005, time constant τNV T = 100 dt = 0.5 τLJ and a PPPMaccuracy of 10−4.To calculate the pressure and energy profiles, the simulation box was dividedinto Nbins parts along the non-periodic z-direction. For every particle in a bin, anindividual per-atom virial tensor σ iab was calculated and averaged across each bin,where a,b could be any of x,y,z and i is the index of the particle. The per-atomvirial tensor in the case of linear polymers with long-range Coulombic interactionsis defined as:σab = − [mvavb+ 12Np∑n=1(r1aF1b + r2aF2b)++12Nb∑n=1(r1aF1b + r2aF2b)+Kspace(ria ,Fib)](3.12)2One way to improve the given model would be to use a hard wall, with polarizable charge onthe electrode to model the polarization of the electrode and formation of the double layer wouldbe the best option. It could be done a methodology described in https://github.com/zhenxingwang/lammps-conp. Without polarization of the electrode, results obtianed via a hard wall should becomparable with the reflective boundary conditions53The first term is a kinetic energy contribution for the given atom. The second termis a pairwise energy contribution, r1 and r2 are the positions of the two atoms inthe pairwise interaction where Np is the total number of neighbors, and F1 and F2are the forces on the two atoms resulting from the pairwise interaction. The thirdterm is a contribution from the covalent bonds, and the last term accounts for long-range Coulombic interactions. This tensor has units of energy and described indetail by Thompson et al. [73]. Two sets of pressure profiles were calculated dur-ing the simulation run. A lateral pressure profile via Plat(z) = −(σ ixx +σ iyy)/2V z jbinand the pressure component along the z axis Pzz(z) = −(σ izz)/V z jbin. During theproduction run, the pressure profile was recorded every 500 steps and averaged inpost-processing. The energy and concentration profiles were calculated in a similarfashion.To perform the counterion condensation analysis we measured the distancesbetween the gel monomers and counterions. If counterions were within rcut = 2σof the gel backbone then it was considered to be condensed. The value of rcutwas chosen to include the first two peaks of the radial distribution function of thecondensed ions. Counterions that were within rcut of more than one backbonemonomer were counted only once. The results for the fractions of the condensedatoms were time averaged across the simulation.To calculate the electrostatic potential profiles, snapshots were saved every1000 steps and imported into VMD’s built-in PME Poisson solver [3]. The po-tentials of each snapshot were averaged to obtain a potential profile φ(z). Thesecurves were fitted toφ(z) = φ0+A · tanh(−x/x0) (3.13)via the parameters x0, φ0, A. Equation 3.13 arises from solving the Poisson equa-tion in the geometry of two semi-infinite slabs with a potential at the interface(ϕ(z−) = ϕ(z+)).MD simulation of PE diodeOur simulation setup comprises of a box of 20x3x3 unit cells forming a cubic lat-tice of cross-linking nodes. The unit cells have lattice constant l = α (presented inTable 3.2) and are connected by polyelectrolyte chains with Nm = 100 monomers.54A fraction f of monomers is then ionized. The ionized monomers on the left sideof the box bear a negative charge q = −1, whereas the chains on the right side ofthe box bear a positive charge q = +1. After the network is created, Nci coun-terions bearing positive and negative charges are inserted into the left and rightsides of the box, respectively to achieve the global electroneutrality. Monomersin the network are tethered to static nodes at the left and right boundaries of thesimulation box. Periodic boundary conditions (PBC) are employed in lateral di-rections only. Reflective walls are put at the boundaries to prevent atoms fromescaping the box. The box parameters are presented in Table 3.2, and for the caseof f = 0.1, lb = 1 σ are z=−lz =−(10+1) ·37=−407 σ and z= lz = 407 σ . Thesystem is then equilibrated to achieve zero lateral pressure. To do that we subjectthe lateral directions of the system (x and y) to an equilibration protocol in the NPTensemble with zero external pressure with time constant τNPT = 100 dt = 0.5 τLJ .This procedure was repeated for two ionizations f = 0.1 and f = 0.25 (i.e. 10 and25% of the monomers are ionized respectively) and two Bjerrum lengths lb = 1as well as lb = 5σ . The full list of parameters including equilibrium volume Veq,backbone ion concentration c f = f Nmons/Veq, effective backbone ion concentrationce f f = fe f f Nmons/Veq are summarized in Table 3.2.Table 3.2: A summary of the simulation parameters used in the present work.f lb, σ fe f f d, σ α, σ lz, σ Veq, 107 ·σ3 c f , 10−4 ·σ−3 k f , σ−1 k−10 , σ ce f f , 10−4 ·σ−3 ke f f , σ−10.1 1 0.100 5 41.3 826.41 1.35 3.76 0.07 29.24 3.76 0.070.1 5 0.073 80 36.7 734.56 0.92 5.51 0.19 11.04 4.32 0.160.25 1 0.250 40 60.3 1205.44 3.57 3.70 0.07 29.81 3.70 0.070.25 5 0.108 120 42.9 858.17 1.63 8.09 0.23 9.15 3.61 0.15The electrostatic cutoff distance was set to relc = 13 σ , above which longer-range electrostatic interactions are calculated via a Particle-Particle-Particle Mesh(PPPM) Ewald solver for slab systems [52]. All simulations were performed in theNV T ensemble with time constant τNV T = 100 dt = 0.5 τLJ and a PPPM accuracyof 10−3.To calculate the electrostatic field profiles, the simulation box was divided intoNbins parts along the non-periodic z-direction. For every particle in a bin, an in-dividual per-atom Coulombic force F icoul was calculated and averaged across each55bin, where i is the index of the particle. To calculate the Coulombic componentof the force, first a production run was performed and positions of every particlewere saved every 100 timesteps. Next, a post-processing simulation rerun wasperformed, where forces and energies of the particles were calculated. Here, theparameters of the LJ and FENE interactions were set to zero in order to isolate theCoulombic component. This is a commonly used trick to calculate the electrostaticfield profiles, that provides a possibility to obtain the field without integrating thecharge density profiles explicitly.We note here that the computational models used to model the sensor and thediode are based on similar models used by [14, 16]. To further verify the validityof the model we performed a benchmark test of the model with a simple electrolyteplaced in a box described in Section A.1.3.2 Continuum theory and modellingTo solve the Poisson-Boltzmann equation (Equation 2.12), a finite-difference methodis used. The domain of interest is split into two and due to the asymmetry of theproblem, only the right half-domain (0≤ z≤ L/2) is taken into consideration anddiscretized into N = 1000 points. Mixed boundary conditions are employed onthe half-domain boundaries. Following Yamamoto and Doi [77], Dirichlet bound-ary conditions are employed on the left side of the half-domain (z = 0) so theelectrical potential is set to zero (ϕ(z = 0) = 0) at the interface. On the rightside of the half-domain (z = L/2), the electrical field is set to a specified valueE = Eboundary. An open source distribution of a PB solver is used which is adaptedfrom https://github.com/pv/scikits.bvp1lg. The module results were first checkedagainst MD simulations as well as analytical solutions for the case of a simpleelectrolyte between electrodes. The MD model, as well as the PB solver has beenbenchmarked on a simple electrolyte system described in Section A.1.56Figure 3.6: Schematic representation of the p-type gel synthesis.Figure 3.7: Schematic representation of the gel synthesis.3.3 Experimental methods3.3.1 Fabrication of PGDPolyelectrolyte gel synthesisTo synthesize3 a polyelectrolyte gel network, we first modified our base materialcrystalline nanocellulose (CNC) to obtain desired level of surface ionization, fol-lowed by a cross-linking of chains into a network. Below we will describe thesynthesis procedure behind creation of an p-type and PE chains, and their physicalcross-linking into a network.In experimental literature, a network with negatively charged PE chains andpositive floating counterions is referred to as a p-type gel. Cellulose nanocrystalswere modified by periodate oxidation and then treated with sodium bi-sulfite toyield the corresponding C2,3 sulfonates as it is shown schematically in Figure 3.6.43This work is done with two UBC graduate students Parya Keyvani (UBC, Department of Chem-ical and Biological Engineering) and Kudzi Nyamayaro (UBC, Department of Chemistry)4Cellulose nanocrystals (CNCs) were purchased from CelluForce Inc. (Montreal, Quebec,Canada) in the form of spray-dried sodium salts with a sulfur content of 0.89 wt.% fraction as re-ported by the supplier.57Figure 3.8: Schematic representation of the electrode setup.A network with positively charged PE chains (also known as poly-cationic)and negative floating counterions are referred to as gel. The nanocrystalline cel-lulose crystals cationized through the reaction between CNCs prepared by sulfu-ric acid and epoxy-propyl-trimethyl-ammonium (EPTMAC). EPTMAC was em-ployed as the cationization agent, yielding cellulose nanocrystals with cationichydroxypropyltrimethylammoniumchloride(HPTMAC-CNC) substituents on the sur-face as shown schematically in Figure 3.7. The aqueous suspension of CNC wasused in the reaction in order to prevent the change in morphology and dimension.5To cross-link the chains into a network we first prepared p-type and CNC so-lutions using de-ionized water under continuous stirring. A weight fraction of 2.5,5 and 7 % of CNC was used. After that, a cross-linking agent agarose with 3weight % was added to them and then the mixtures were boiled on a hot plate un-der continuous stirring. The mixtures were cooled down to room temperature. As aresult the chains in the solution gelated which was verified via relevant rheologicalmeasurements.The charged groups, on the surface of the material were quantified using con-ductometric titration and elemental analysis; results are summarised in Table 3.3.The charge density determined by conductometric titration for pCNC was found to5Poly (vinyl alcohol) MW= 146,000-186,000, 99% hydrolyzed ,agarose type I, Glycidyltrimethylammonium chloride (2,3-Epoxypropyl)trimethylammonium chloride)58be 1.4010 mmol g−1 and was significantly higher than CNC at 0.1986 mmol g−1.The same trend was observed in the sulphur content determined by elemental anal-ysis, which verifies the increase in sulfonate groups [81]. The surface charge ofnCNC was determined to be 0.1305 mmol g−1 and was confirmed by the nitrogencontent elemental analysis. In all cases the charge density values from conducto-metric titration are lower than EA results, possibly due to the underestimation ofthe charge content associated with the use of mixed bed resin .Elemental analysis (EA) was performed at the UBC Mass Spectrometry andMicroanalytical Laboratory on a Fisons Instruments elemental analyzer EA 1108using flash combustion. The surface charge was determined by conductometrictitration, by employing previously reported methods. Briefly, approximately 0.1g of dried nCNC was suspended in 50 ml of 1 mM NaCl in a 100 ml three-neckflask. The conductivity of the suspension was monitored during the titration with10.11 mM AgCl from a dropper set to disperse 0.1 mL/min. The sulphate contentin pCNC was determined by titrating 0.1 g of pCNC in 50 ml of 1 mM NaCl with9.980 mM NaOH.Rheology was used to evaluate the strength and stability of the agarose gelsdoped with the CNC, nCNC and pCNC. Strain sweep oscillatory shear tests high-lighted that the linear viscoelastic response of the gels remained relatively constant.However, the critical strain for the filled hydrogels were slightly lower which indi-cates that the doped hydrogel is more fragile than the pure agarose hydrogel. Thishas been shown to be the case for high loadings of CNC; the particles disrupt theagarose network, reducing the connectivity and hence lowering elasticity.Sample ζ CI Titration charge density (mmol g−1) EA charge density (mmol g−1)CNC -47.7 ± 2.11 81 0.1986 ± 0.0033 0.2047 ± 0.0109nCNC +49.8 ± 0.89 60 0.1305 ± 0.0004 0.4607 ± 0.0036pCNC -38.8 ± 9.66 45 1.4010 ± 0.0110 1.6328 ± 0.0703Table 3.3: Zeta potential, crystallinity index, charge density from conducto-metric titration and charge density from elemental analysis of unmodifiedCNC, nCNC and pCNC. Error bars represent a standart deviation fromthe mean after 3 measurements59Figure 3.9: Zeta potential from conductometric titration of unmodified CNC,nCNC and pCNC. Error bars represent a standard deviation from themean after 3 measurements.Figure 3.10: Schematic representation of the fabrication method to fabricatea PGD.Assembly of the diode and experimental studyCylindrical slabs were cut from CNC-AG hydrogels and inserted into the 10 mmholes in a non-conductive/non-reactive spacers with thickness of 1.2 mm. Thespacers with agarose gels containing either pCNC or nCNC were put into con-tact and then clamped between two ITO electrodes. The characteristic I–V curveswere obtained using a computer controlled Solarton 1287A potentiostat/galvano-stat, Hampshire, UK. The performance of the agarose CNC based diode was in-vestigated in terms of type of applied bias (step or sweep shown in Figure 3.11),sweep scan rate, gel thickness and nCNC/pCNC concentration. All procedures60Figure 3.11: Schematic representation of study protocol of a PGD. A lin-ear sweep between voltage range at different rate(black), and a stepbias(red).were carried out at room temperature (295 ± 2 K).61Chapter 4Piezoionic effectOur goal for studying the piezo-ionic effect using Molecular Dynamics was to an-alyze the connection between the system’s inhomogenuity in ionization and gener-ated potential and pressure difference. So, we begin our investigation of the piezo-ionic effect by analyzing a system with ionization degrees f1 = 0.1 and f2 = 1.0,i.e., 10% of monomers are ionized on the left and 100% on the right side of thematerial. We run our simulations at different Bjerrum lengths lb achieved by vary-ing the temperature of the system T . Thus we first obtain an qualitative picture bystudying the snapshots of the system, followed by detailed analysis using Manningcounterion condensation criterion and studying the potential energy differences be-tween two parts of the device. We finish our investigation by investigating the con-nection between ionization difference and the generated Nernst-Donnan potentialsas well as the pressure gradient between two parts.4.1 Three regimes of Donnan equilibriumIn Figure 4.1 we show the conformation of the gel and counterions, by lookingat several snapshots of the system at three different temperatures T = 1.1,0.3,0.1(lB = 0.9,3.3,10.0). A qualitative difference is observed in the behavior at high,medium and low temperatures. At high temperatures (regime 1 or weak-weakregime), both sides of the gel are in the limit of weak coupling where entropicforces dominate the interactions. Here counterions float uniformly within the net-62Figure 4.1: Snapshots of a system with ionization degrees f1 = 0.1 andf2 = 1.0 at three temperatures T = 1.1,0.3,0.1 (at lB = 0.9,3.3,10).In snapshots white balls are: uncharged backbone monomers, red: neg-atively charged monomers, blue: tethered static monomers, green: pos-itively charged counterions. For illustration purposes, only counterionswithin the dashed box are shown. Regime 1 (weak-weak): weak cou-pling for both sides of the sensor with no counterion condensation onboth sides Regime 2 (weak-strong): weak coupling on the left side, andstrong coupling on the right side of the sensor Regime 3 (strong-strong):strong coupling regime for both parts of the systemwork. They exhibit high osmotic pressure and cause significant swelling. Due tothe high degree of swelling, the boundary between two parts of the gel is pushedfurther left. As the temperature decreases to T = 0.3 (lB = 3.3 σ ) the system tran-sitions into regime 2 (weak-strong). The particles on the left side of the box remainin the limit of weak coupling and continue to behave as they did in regime 1, whilethe counterions on the right side undergo counterion condensation and transition to63the strong-coupling regime. Positively charged floating counterions become boundto the negatively charged backbone macroions. The bound counterions no longercontribute to the osmotic pressure, hence reducing the levels of swelling of thenetwork. When the temperature becomes even lower (regime 3 or strong-strong),counterions in both parts of the system transition to the strong coupling regime.They condense and become bound to the backbone macroions. Since the fractionof the free-floating counterions is further reduced due to the condensation on bothparts, the exerted osmotic pressure decreases even further. This causes significantlylower degrees of swelling than in the case of regime 1, as seen in Figure 4.1.These trends can be rationalized with Manning criterion for counterion conden-sation described in Section 2.1.3. The theory compares the inter-charge spacing inthe backbone of the PE to the Bjerrum length lb of the system. When the spacing abecomes smaller than the lb counterions condense. For the present gel, the spacingof charges in gel 1 is a1 = Ree/Nm f1 ' 5, while a2 ' 1. In regime 1, the Bjer-rum length lB = 1/T = 0.91 is smaller than both a1 and a2, hence no counterionscondense. In regime 2, however, lB = 3.33 exceeds a2 but not a1, so counterionscondense in gel 3 only. Finally in regime 3, lB = 10 exceeds both a1 and a2 andcondensation occurs in both gels.To quantify the counterion condensation, we performed a counterion conden-sation analysis and compared it to predictions of Manning criterion described inSection 2.1.3. In Figure 4.2, we show the fraction of free counterions ffree,i =fi(1−Qi), where Qi is the fraction of condensed ions. The theory predicts that athigh temperatures no counterion condensation will occur and the effective ioniza-tion ffree,i should remain equal to the actual ionization fi. As the temperature islowered, counterions are expected to condense, hence ffree,i should decrease. Asexpected, we observe that the fraction of free counterions saturates at high tem-peratures, indicating that the electrostatic forces are not strong enough to facilitatecondensation. Also, at a certain temperature (which depends on the ionization ofthe gel), the gel undergoes a transition where the effective ionization starts to de-crease with decreasing temperature. The three regimes of the Donnan equilibriumintroduced previously in Figure 4.1 can also be observed depending on the onsetof counterion condensation in each part of the simulation box. As the temperaturedecreases from T=1.1, where ions in both parts of the gel are not condensed, to640.00.20.40.6f free,10.00.51.0f free,20.0 0.2 0.4 0.6 0.8 1.0 1.2T [ ]0.51.0f free,2f free,1Regime 3 Regime 2 Regime 1Figure 4.2: Analysis of the fraction of free counterions (open symbols) com-pared to the Manning theory (filled symbols) for the gel slabs with ion-ization degrees f1 = 0.1 and f2 = 1.0 (blue) and f1 = 0.25 and f2 = 0.5(cyan). A: Fraction of free counterions f f ree,1 = f1(1−Q1) on the leftside of the gel. B: Fraction of free counterions f f ree,2 = f2(1−Q2) onthe right side of the gel. C: Ratio f f ree,2/ f f ree,1. Vertical dashed linescorrespond to the temperatures where counterion condensation occursin one part of the gel as indicated by the change in the curve behaviour.T=0.7 the right side of the gel undergoes counterion condensation as indicated bythe decreasing fraction of counterions, see Figure 4.2 B and the system transitionsfrom regime 1 to regime 2. When the temperature is further decreased to T=0.5,the counterions on the left side of the gel also condense, as indicated by the change65in curve behavior in Figure 4.2 A. The three regimes can also be observed on thegraph of the ratio of effective ionization degrees f f ree,2/ f f ree,2, see Figure 4.2 C.In Figure 4.2, we also compare this effective ionization ffree,i that was calcu-lated by directly analyzing the number of condensed counterions with the Manningprediction given by Equation 2.9. The simulation curves are in a excellent agree-ment with this simple prediction in case of low ionization, however, the discrepan-cies increase in the case of higher ionizations (bottom panel). The discrepanciesare possibly due to the fact that the Manning model was derived from solving aPB equation for a straight rod. However, in case of high ionization the system hasnumerous charged particles around the nodes, which have a star-topology ratherthan one of a straight chain.We note here that altering lb in experimental setting is rather cumbersome.Change in temperature often leads to counter-balancing change in dielectric permi-tivity of the medium rendering the Bjerrum length lb ∼ 1/εrT to remain constant.Experimentally, counterion condensation is induced not via changing the Bjerrumlength, but via changing the ionization of the system and keeping the lb constant.For example, one could use the a poly-(AA-co-AM) copolymer gels and changethe fraction of AM (a neutral monomer). When fraction of AA (a polyelectrolyte)becomes large enough, distance between the units becomes small enough (smallerthan lb) the resulting system will collapse [54].4.2 Potential energiesThe physical transition from regimes 1 to 3 shown in Figure 4.1 also alters thepotential energy profiles shown in Figure 4.3. Here we plot the normalized profilesof the pair energy of interaction along the z-axis. We split the pair energy intoLennard-Jones (repulsion due to the excluded volume interaction) and Coulombiccomponents. To be able to observe the qualitative changes in the energy of thesystem with temperature, the energy values were normalized by temperature. Theenergy profiles in Figure 4.3 exhibit qualitatively different behavior across the threeaforementioned regimes. In regime 1 (high temperatures T = 1.1− 0.5), both theleft and the right part of the normalized energy profiles show very little change withtemperature and approach a common curve. In this regime, the energy is essentially66proportional to temperature. In regime 2 (medium temperatures T = 0.4−0.2), theleft part of the system shows no change with temperature, but the right one exhibitsan increase in the value of the Lennard-Jones-component and a decrease in theCoulombic component. Finally, in regime 3 (T = 0.1), the energy profiles in bothparts change. Instead of a evenly distributed gas of counterions, all counterions arenow concentrated in the vicinity of the backbone ions. This drastically increasesthe repulsive Lennard-Jones component and decreases the Coulombic componentof the energy.0.2 0.4 0.6 0.850Eci coul/(NT)T=0.1T=0.2T=0.3T=0.4T=0.5T=0.7T=0.9T=1.10.2 0.4 0.6 0.8z/lz0.00.1Eci lj/(NT)Figure 4.3: Normalized energy profiles (A: Coulombic contribution, B:Lennard-Jones contribution) of the system with f1 = 0.1, f2 = 1.0 fordifferent regimes. Here N is the total number of atoms in the system.4.3 Nernst-Donnan model for the electrostatic potentialWe start by testing whether the Nernst-Donnan picture describes the molecularsimulations correctly. Nernst-Donnan equation connects ionization difference topotential difference observed between parts of the gel, so we commence by com-puting the ion concentration profiles in Figure 4.4. With increasing temperature,the counterion profile becomes smoother as the itinerant ions equilibrate more eas-ily across the system against the electrostatic energy penalty. We denote the con-centrations in the two gels as c1 and c2, respectively. In order to systematicallydefine these values in the two parts of the gel, we take the peak value of the con-670.2 0.4 0.6 0.8z/lz0.51.01.52.02.5cci[1033 ]T=0.1T=0.2T=0.3T=0.5T=0.7T=0.9T=1.1Figure 4.4: Counterion concentration profiles of the system with f1 = 0.1,f2 = 1.0 for different temperature regimes. Colors are the same as inFigure 4.3centrations away from any boundaries. Just as before, we introduce a value cfree,iwhich is a concentration of uncondensed counterions cfree,i = ci(1−Qi). Equa-tion 2.11 predicts that the potential scales linearly with temperature with the coef-ficient of the proportionality given by ln(cfree,2/cfree,1). To test this, the potentialenergies q∆ϕ of four different gels are scaled by ln(cfree,2/cfree,1) and plotted as afunction of temperature, see Figure 4.5. All data sets collapse reasonably well ontoa linear master curve with slope unity (circles).Additionally, the master curve has zero offset along the vertical axis, with aslight divergence from linear behaviour at low temperatures. This confirms thatexcluded volume interactions and hence the excess chemical potential in Equa-tion 2.11 are negligible in the current model, especially in the high temperatureregime.In experiments, the degree of ionization is not directly controlled, but is as-sumed to respond to external pressure changes. The electric potential can be con-nected to the applied pressure in two steps: understanding the behavior of the pres-sure profiles and then connecting the potential to the observed pressure differencebetween the two gels. To understand the behaviour of the pressure profiles at dif-ferent temperatures, the pressure components Plat and Pzz horizontal and vertical tothe gel slab are computed as a function of the z-coordinate. The lateral pressure Platshows a step-like behavior with higher pressure on the higher ionized side, while680.0 0.2 0.4 0.6 0.8 1.0T [ ]0.00.20.40.60.81.01.2q*[]Figure 4.5: Normalized electrostatic potential as a function of temperature.Circles: q∆ϕ∗ = q∆ϕ/ ln(cfree,2/cfree,1). Red: f1 = 0.1, f2 = 0.2; ma-genta: f1 = 0.1, f2 = 0.3; green: f1 = 0.1, f2 = 0.5; blue: f1 = 0.1,f2 = 1.0. The dashed line has slope unity.Pzz remains constant across the gels as required by the condition of mechanicalequilibrium. In Figure 4.6, we plot the potential difference between the gel partsq∆ϕ = q(ϕ2−ϕ1) as a function of lateral pressure difference Plat,2−Plat,1. Thefigure shows that the data points fall nearly perfectly onto a straight line, especiallyat lower temperatures. The gels with higher ionization difference f2− f1 exhibithigher pressure difference values and consequently higher potential values. Thisresult shows that the energy of expansion scales linearly with the electrostatic en-ergy change, which can also be inferred from the equality of chemical potentialsthat reads ∆P = ∆Π. Specifically, the minimum property of the Gibbs free energythat reads V∆P= N∆µ , leads to ∆P∼ ∆Π∼ cfree,1kBT ln(cfree,2/cfree,1), which im-plies a linear scaling between pressure and Nernst potential when combined withq∆φ = kBT ln(cfree,2/cfree,1) from Figure 4.5.690.5 1.0 1.5 2.0 2.5 3.0(P2 P1)/cfree, 1[ ]0.00.51.01.5q[]Figure 4.6: Potential difference as a function of the lateral pressure differ-ence. Red: f1 = 0.1, f2 = 0.2; magenta: f1 = 0.1, f2 = 0.3; green:f1 = 0.1, f2 = 0.5; blue: f1 = 0.1, f2 = 1.0. Here V is the volume ofthe simulation box.4.4 Summary• We have identified a linear relationship between the effective concentrationdifference (effective ionization) and the Donnan potentialq∆ϕ∗ = q∆ϕ/ ln(cfree,2/cfree,1) = T• The model remains valid in the limit of weak coupling (i.e., high Bjerrumlengths) if accounted for counterion condensation• Unlike the experimental setting where applied pressure caused counteriondissociation and buildup of Donnan potential, an MD model where ioniza-tion difference is setup by design results in a buildup of electrostatic poten-tial, as well as a pressure gradient between two parts of the gel• These results provide a framework for the analysis of experimental measure-ments of PE soft sensors in the limit of low salt concentrations (large Debyelengths ).70Chapter 5Molecular Dynamics modeling ofPE diodesTo create a reliable PB model for PGDs that is consistent with the MD results, wefirst need to understand the behavior of PE gels at different ionizations f , differentBjerrum lengths lb and strengths of the applied applied electric field E. This will befollowed by a comparison between direct MD results and the PB model in differentcoupling regimes. We note again that the present study focuses on the electrostaticsin model PGDs without explicit water molecules, does not consider electrolysisat the boundaries of the diode, i.e., the regime with suppressed currents at theelectrodes.5.1 Weakly ionized chain with f = 0.15.1.1 Weak coupling (lb = 1)We start our investigation by demonstrating that for a weakly ionized chain f = 0.1in the limit of weak coupling (lb = 1) a PB model provides a good description of theMD results. We illustrate that by comparing the charge density profiles obtainedfrom the MD simulations and the PB model, followed by the analysis of the electricfield profiles.In Figure 5.1 we show the spatial charge density profiles at different strengths71Backbone ions“+” counter ions“-” counter ionsPB modelMD simulationEboundary = 0Eboundary = -0.1Eboundary = 0.1f = 0.1, lb=1Figure 5.1: Ion spatial density profile of the system with f = 0.1, lb = 1 atdifferent strengths of applied electric field E = 0 (top), E =−0.1 (mid-dle), E = 0.1 (bottom). Blue: negatively charged itinerant ions; red:positively charged itinerant ions; green: gel backbone ions (both pos-itively and negatively charged); solid line: PB model solution; dashedline: MD profiles. Error bars represent ±1σstd confidence intervals (forclarity purposes only every 6th error bar is shown).of the applied field. Dashed lines in Figure 5.1 correspond to MD data and thesolid lines represent the solution of the PB model. The blue and red colors corre-spond to itinerant ions bearing negative or positive charge, respectively, and greencorresponds to the backbone ions. The top panel of the Figure 5.1 represents thescenario when no electric field is applied E = 0, the middle panel represents theblocking regime E = −0.1, and the bottom panel represents the rectifying regimewith E = 0.1. We see that regardless of the electric field applied the backboneion concentration is considered to be constant c f , which is a function of the diode72PB modelMD simEboundary = 0Eboundary = -0.1Eboundary = 0.1f = 0.1,lb=1Figure 5.2: Spatial profile of the electric field inside of the diode at differ-ent strengths of applied field (black: E = 0, blue: E = −0.1, red:E =+0.1) for the weakly ionized system in the limit of weak coupling( f = 0.1, lb = 1). Dashed lines represent the MD data; solid lines arethe corresponding PB model curves. Error bars represent ±1σstd confi-dence intervals (for clarity purposes only every 6th error bar is shown).ionization and monomer density. Unlike the backbone ions, the itinerant (float-ing) ions are free to move throughout the gel. When no electric field is appliedthe counterion profile is constant and follows the backbone ion profile on one sideof the diode c f and is equal to 0 on the other. At the interface between the twosides of the diode (z= 0), the charge density decays following Boltzmann distribu-tion creating a depletion zone 2 · lGC thick, where lGC = k−10 = (8pilbc0)−1/2 is theGouy-Chapman length as shown in Figure 1.6. When a negative field is applied nochanges are observed in the bulk and at the interface of the diode, but an increaseis observed in the counterion concentration on the boundaries (z/lz = ±1). Ap-plying a positive electric field has an opposite effect as it decreases the counterionconcentration on the boundaries.Regardless of the electric field applied, MD data shows a periodic structurewith a period α corresponding to the spacing of the cubic lattice of cross-linkingnodes.Thee large oscillations are an artifact of the 6-fold symmetry of the cubiclattice. Cross-linking in real gels is usually more random and hence has less sym-73metry. Notwithstanding the oscillations, the present results are relevant as they pro-vide a good model for studying the interface between two systems and were alsoused other similar simulations [14]. Unlike the backbone ions that are constrainedin a network via bonds, the itinerant (floating) ions are free to move throughoutthe gel. This results in much smoother counterion profiles. The concentration fluc-tuates around an average value c f = f · c in the bulk of the diode, but decays atthe interface following a Boltzmann distribution to the value c0 just as in the PBmodel. Despite being constrained in a network, the backbone ions still have theflexibility to interact and move, although to a lesser degree than the itinerant ions.Due to this interaction, the gels at the interface attract each other and coalesce,creating effectively two interfaces in the middle of the diode as can be seen in Fig-ure 5.3. The coalescence pulls the gels from the boundaries inward and decreasesthe concentration of backbone and itinerant ions at (z/lz = ±1). The coalescenceof the gels also results in a peak of the backbone density profile in the counteriondepletion zone at the interface as shown in Figure 5.1.In the case of zero external field (Figure 5.1 top panel), the results from thePB model and MD simulations generally agree throughout the diode except at theboundaries and at the interface. In the bulk region of the diode, both itinerant andbackbone ions show similar behaviour in the MD and PB models. At the interfacebetween the gels (center of the diode), the MD density profiles for itinerant ionsmatch the PB model, however the MD backbone ion density values are signifi-cantly higher than the PB prediction due to the gel coalescence described aboveand illustrated in Figure 5.3. On the outside boundaries (z/lz = ±1) of the diodewe see constant density values for the PB model, and a decrease in density for theMD simulation.When a negative electric field is applied (blocking regime of the PGD), the MDmodel shows that the ions are behaving differently at the boundaries than the PBmodel would predict (Figure 5.1 middle panel). On the left side of the diode thePB model shows a significant increase in concentration of positively charged ionscompared to only a slight increase in the MD model, while the concentration of thenegatively charged itinerant ions are zero for both PB and MD model.The agreement between the PB and MD results becomes better when the sys-tem is subjected to a positive external field (rectifying regime, Figure 5.1 bottom74Figure 5.3: Ion spatial density profile of the system with f = 0.25, lb = 1 atdifferent strengths of applied electric field E = 0 (top), E =−0.1 (mid-dle), E = 0.1 (bottom). Blue: negatively charged itinerant ions; red:positively charged itinerant ions; green: gel backbone ions (both pos-itively and negatively charged); solid line: PB model solution; dashedline: MD profiles. Error bars represent ±1σstd confidence intervals (forclarity purposes only every 6th error bar is shown).panel). On the left side of the diode both the PB model and the MD simulationproduce a decrease in the concentration of the positive ions and an increase in theconcentration of the negative ions that have migrated from the opposite side ofthe diode. A similar picture is observed on the right side of the diode due to theasymmetry of the system.Surprisingly, none of the discrepancies in the concentration profiles of the MDand PB models result in significant differences in the electric field profiles as shownin Figure 5.2. The figure presents the spatial profile of the electric field inside of the75diode at different strengths of applied field (shown in different colors) with dashedlines corresponding to the MD data and the solid lines to the PB model curves.When no external field is applied, both the PB and MD curves show zero electricfield at the diode boundaries and in the bulk of the gel, and have a negative dip atthe interface between the two gels. However, the PB model fails to reproduce thetriple peak structure at the diode interface which occurs due to the collapse of thegel network at the interface. When an external field is applied, both PB and MDmodels show that the electric field curves remain similar to the no-field regimethroughout the gel, except at the boundaries where it experiences an increase ordecrease depending on the sign of the applied field.The comparison of PB and MD results provides great insight into the behaviourof the diode on its boundaries, in the bulk as well as at the interface. On the bound-aries, PB and MD models result in similar electric field profiles, despite havingdifferent ion concentration. The reason behind this mismatch, perhaps, can be seenby considering the details of the PB and MD models. In the case of the PB model,the electric field at the boundary has to be equal to the specified value as per theVon-Neumann boundary condition. Since the backbone ion concentration is set toa constant (c = c f ) value throughout the gel, the concentration of the itinerant ionsis driven solely by the sign of the external field and the value of the Gouy-Chapmanlength k−10 .In the case of the MD simulation, however, the backbone ion concentrationis not constant but fluctuates and is able to change due to the flexibility of thenetwork. Additionally, the electric field in the MD simulation appears explicitly asan additional force on the ions. The network flexibility, as well as application ofthe external field produce a system where concentration profiles are driven firstlyvia the change in the backbone concentration, and only then by the counterionequilibration (determined by the sign of the field and the Gouy-Chapman length).The counterion equilibration happens in a manner so they screen the electric fieldin the bulk of the gel (0 |z|  lz) indicating ”conductor-like” behaviour andproduce a field Eboundary at the interface. For example, for a negative externalelectric field, counterions are pushed out of the gel to the box boundary. Sincethe backbone density in the case of MD is lower than the backbone density in thePB model, the counterion density at the boundary is also smaller. In the case of76positive external electric field, itinerant ions are pushed inside the gel. Therefore,despite having different backbone ion profiles counterions equilibrate themselvesin a way that produces curves similar to the PB model.Finally, the most interesting part of the profile is the interface between thetwo gels. Here both the PB and MD models produce a negative dip of the electricfield. Unlike the PB model, however, which produced a single peak, the MD modelshows a triple peak structure. A single peak arises from the incomplete chargeneutralization in the transition region between the gels, where counterions spillfrom their respective gels into the interface regions and recombine with counterionsof the opposite charge as shown in the sketch presented in Figure 1.6. The widthof the peak is determined by the Gouy-Chapman length k−10 = (8pilbc0)−1/2. Thetriple peak in the electric field profile comes from the gel chain coalescence atthe interface. Unlike the case of a constant backbone concentration (used in thePB model), that has a single interface, the coalescing gels create two interfacesbetween the narrow collapsed region and the bulk of the two gels as seen from thetop panel of Figure 5.2.5.1.2 Strong coupling (lb = 5)Having understood the behaviour of the system in the weak coupling regime, wenow increase the Bjerrum length and consider the regime of strong coupling. Asdiscussed before, in this regime a fraction of itinerant ions becomes strongly cou-pled to their respective backbone ions, reducing effective ionization of the gel (de-scribed in Equation 2.9 and summarized in Table 3.2) and producing different con-centration and electric field profiles. For all values of the electric field examinedin Figure 5.4, the PB curves now exhibit a narrower and steeper interface as theGouy-Chapman length has decreased. As a result of charge conservation, we alsosee smaller changes at the boundaries when an external field is applied.The MD simulation data shown as dashed lines presents a different picture.On the gel boundaries, at zero electric field we see a scenario similar to the caseof weak coupling. When a positive or a negative bias is applied, we see that thegel concentration on the boundaries changes. It decreases when the negative fieldis applied, and increases when a positive field is applied. The counterions follow77the concentration of the backbone gel ions. The concentration profiles in the bulkof the diode differ from the case of weak coupling as they show reduced oscilla-tions closer to the diode boundaries, suggesting that the gel is being ”pulled” awayfrom the center. The profiles for the itinerant ions near peaks (blue and red) practi-cally mirror the concentration of the backbone ions (green), confirming counterioncondensation. We also notice an increased bulk concentration of the gel c f , sug-gesting a decreased volume of the simulation, and hence a decreased equilibriuminter-crosslink distance. The distance is reduced as the equilibrium PE chain sizeis reduced following Ree = Nmb f 1−νe f f < Nmb f1−ν due to counterion condensation.At the interface region between two gels we see a much narrower depletion zone,and higher peaks for the concentration profiles for the itinerant ions. It is worthmentioning, that in the collapsed gel the counter ions condensed on the backboneare shared by two chains at the same time.The more complex behaviour of the diode in the limit of strong coupling resultsin different electric field profiles shown in Figure 5.5. The PB curves in all casesremain similar to the weak coupling regime throughout the gel except near the in-terface in the middle of the diode, where the depth and width become narrowerand deeper (Figure 5.5). The MD profiles, however, are now qualitatively differ-ent than the PB solution. They show no dip behaviour at the interface. The MDsimulation thus predicts that PGDs cease to display rectifying behavior at higherBjerrum lengths, an effect not predicted by the PB treatment.An in-depth comparison between the PB model and the MD simulation shedsmore light into the behaviour of the gel in the limit of strong coupling as well ashighlights the limitations of the PB model. In this regime, the PB and MD resultsare different at the interface. In the PB model the backbone ion concentration at theinterface of the diode changes instantaneously from negative to positive ionizationwith a constant absolute value. As a result, the concentration profiles and electricfields are driven primarily by the itinerant ion equilibration. The MD model ac-counts for many more details not considered in the PB model: the backbone ionconcentration is not static due to the periodicity of the gel as well as field inducedgel deformation. In the limit of strong coupling, counterions in the MD modelcondense and thus reduce the effective ionization of the system. Moreover, theinterface region in the MD model is not abrupt as in the PB model but is gradual78Backbone ions“+” counter ions“-” counter ionsPB modelMD simulationEboundary = 0Eboundary = -0.1Eboundary = 0.1f = 0.1, lb=5Figure 5.4: Ion spatial density profile of the system with f = 0.1, lb = 5 atdifferent strengths of applied electric field E = 0 (top), E =−0.1 (mid-dle), E = 0.1 (bottom). Blue: negatively charged itinerant ions; red:positively charged itinerant ions; green: gel backbone ions (both pos-itively and negatively charged); solid line: PB model solution; dashedline: MD profiles. Error bars represent ±1σstd confidence intervals (forclarity purposes only every 6th error bar is shown).due to a strong collapse of the PE network at the interface shown in Figure 5.3.Due to the effective reduction in the fraction of mobile counterions, and the rampup of the backbone ion concentration the dip in the electric field (and hence thepotential difference between the gels) is now much smaller. Hence, to improve thePB model one does not only need to consider counterion condensation but also amore realistic concentration profile of the backbone ions.79PB modelMD simEboundary = 0Eboundary = -0.1Eboundary = 0.1f = 0.1,lb=5Figure 5.5: Spatial profile of the electric field inside of the diode at differ-ent strengths of applied field (black: E = 0, blue: E = −0.1, red:E =+0.1) of the weakly ionized system in the limit of strong coupling( f = 0.1, lb = 5). Dashed lines represent the MD data; solid lines arethe corresponding PB model curves. Error bars represent ±1σstd confi-dence intervals (for clarity purposes only every 6th error bar is shown).5.2 Effect of ionization ( f = 0.25)So far we have studied the electric field profiles of a system with every 10th of itsmonomers bearing a charge. By varying the Bjerrum length lb we found that theelectric field profiles change when the lb becomes close to the inter-charge distancein the chain. Higher ionizations in principle may also lead to discrepancies betweenPB and MD results even for low Bjerrum lengths. To test that, we have used the sys-tem studied in 5.1 and increased the fraction of charge-bearing monomers to 25 %.Figures 5.6 and 5.7 show the concentration and electric field profiles, respectively,of this system in the limit of weak coupling (lb = 1). As before, the dashed linesrepresent MD results and solid lines represent the PB solution. Overall, the resultsfollow a similar trend as in the case of lower ionization. The concentration profilesfor the regime of weak coupling qualitatively follow the scenario of low ionization,except in the bulk region where the bulk backbone density as well as its oscillationsare higher due to an increased ionization fraction.80Backbone ions“+” counter ions“-” counter ionsPB modelMD simulationEboundary = 0Eboundary = -0.1Eboundary = 0.1f = 0.25, lb=1Figure 5.6: Ion spatial density profile of the system with f = 0.25, lb = 1 atdifferent strengths of applied electric field E = 0 (top), E =−0.1 (mid-dle), E = 0.1 (bottom). Blue: negatively charged itinerant ions; red:positively charged itinerant ions; green: gel backbone ions (both pos-itively and negatively charged); solid line: PB model solution; dashedline: MD profiles. Error bars represent ±1σstd confidence intervals (forclarity purposes only every 6th error bar is shown).The electric field profiles are also similar to the case of low ionization, but onlyat the boundaries and the bulk of the diode. The PB solution for f = 0.25 shows awider and deeper interface due to an increased concentration of charges. The MDresults, however, show an even wider and shallower profile, and the magnitude ofthe dip has been reduced by about a factor of two. The discrepancy originates fromthe attractive forces present between the two parts of the gel as well as their flexi-bility. At higher charge density, the attractive forces between two sides of the diodeare stronger which results in a larger degree of coalescence between two parts. The81PB modelMD simEboundary = 0Eboundary = -0.1Eboundary = 0.1f = 0.25,lb=1Figure 5.7: Spatial profile of the electric field inside of the diode at differentstrengths of applied field (black: E = 0, blue: E =−0.1, red: E =+0.1)in the limit of weak coupling (lb = 1) and ionization fraction f = 0.25.Dashed lines represent the MD data; solid lines are the correspondingPB model curves. Error bars represent ±1σstd confidence intervals (forclarity purposes only every 6th error bar is shown).collapsed chains at the interface produce a large zone where the polyelectrolytechains cancel each other, as well as share their respective counterions. This pro-duces a wider interface between the gels with a gradual charge density decrease.The ramp in the backbone concentration at the interface is fundamentally differentfrom the boundary condition used in the PB model, where the interface is treatedas abrupt.In the limit of strong coupling, chains collapse at the interface becomes evenstronger producing a larger electrically neutral zone, and further reducing the dipin the electric field profile shown in Figure 5.9. The almost-zero dip at the interfacemeans that the system with 25% of its ions ionized in the limit of strong couplingexhibits no rectifying behaviour and acts like a conducting electrolyte with zerofield inside it.82Backbone ions“+” counter ions“-” counter ionsPB modelMD simulationEboundary = 0Eboundary = -0.1Eboundary = 0.1f = 0.25, lb=5Figure 5.8: Ion spatial density profile of the system with f = 0.25, lb = 5 atdifferent strengths of applied electric field E = 0 (top), E =−0.1 (mid-dle), E = 0.1 (bottom). Blue: negatively charged itinerant ions; red:positively charged itinerant ions; green: gel backbone ions (both pos-itively and negatively charged); solid line: PB model solution; dashedline: MD profiles. Error bars represent ±1σstd confidence intervals (forclarity purposes only every 6th error bar is shown).5.3 Modified Poisson Boltzmann modelIn the standard PB model, the degree of ionization does not depend on the Bjerrumlength and hence does not consider the reduction of ionization due to counterioncondensation. The reduction of effective ionization ( f → fe f f ) can be either ob-tained via Manning theory eq. (2.9), or taken directly from the MD simulations aswas done in [75]. By substituting f with fe f f we also replace k f ,c f → k fe f f , c fe f fin eq. (2.13). We found that solving modifying the PB equation by only substitut-83PB modelMD simEboundary = 0Eboundary = -0.1Eboundary = 0.1f = 0.25,lb=5Figure 5.9: Spatial profile of the electric field inside of the diode at differentstrengths of applied field (black: E = 0, blue: E =−0.1, red: E =+0.1)in the limit of weak coupling (lb = 5) and ionization fraction f = 0.25.Dashed lines represent the MD data; solid lines are the correspondingPB model curves. Error bars represent ±1σstd confidence intervals (forclarity purposes only every 6th error bar is shown).ing k f ,c f → k fe f f , c fe f f reduces the depth of the dip at the interface of the diodeby half, but doesn’t account for an increased width of the profile demonstrated bythe MD simulations. One needs to account additionally for the gradual increasein the backbone charge density at the interface. With this in mind, we also relaxthe assumption of fixed gel concentrations c f , and replace the abrupt interface be-tween them with a smoothed interpolation via a tanh(z/d) function that mimics aramp in backbone ion concentration for the diode interface. This model shown inEquation 5.1 will be referred to as modified PB model:d2ϕ(z)dz2= k2o sinh(ϕ(z))− k2e f f tanh(z/d) (5.1)The smoothing function tanh(z/d) has been previously used to describe thedensity profile of a freely-standing polymer film [68].The parameter d was treatedas a fit parameter in the present work as shown in Table 3.2, but is a function of lb, fand could be obtained from a more detailed analysis of the ramp up in the charge84PB (modif.)MDlb = 5lb = 1f = 0.1f = 0.25Figure 5.10: Spatial profile of the electric field inside of the diode at differ-ent ionizations ( f = 0.1 top, f = 0.25 bottom) and couplings (black:lb = 5, red: lb = 1). Circles represent the MD data; solid lines are thecorresponding modified PB model curves with parameter d describedin Table 3.2. Error bars represent±1σstd confidence intervals (for clar-ity purposes only every 6th error bar is shown).density profile. An illustration of the effect of these modifications is shown in Fig-ure 5.10. In the limit of weak coupling (lb = 1) we see that the modified PB modelaccounts for the small value of d which is the length of the ”transition” regionat the interface and produces narrow and deep dip. As the ionization or Bjerrumlength increases the model produces a shallower and wider interface accounting fora larger transition region accounted for by a larger value of d. For the given valuesof lb,c f ,d the modified PB model can produce a better description of the dip at theinterface of the PGD and hence better predicts its rectifying behaviour (or the lackof thereof).5.4 Summary• A standard PB model is able to predict the ion density and electric field pro-files of the PGD in the regime where currents are surpressed at the electrodes• The model fails to reproduce the regime of strong coupling (i.e., large Bjer-85rum lengths) due to counterion condensation as well as a collapse of polymerchains at the interface• A modified PB model has been proposed that accounts for counterion con-densation (via an effective inverse length scale ke f f ) as well as the chaincollapse (via introduction of the ramp-up parameter tanh(z/d))d2ϕ(z)dz2= k2o sinh(ϕ(z))− k2e f f tanh(z/d)• These results allow to modify the Yamamoto model for better prediction ofthe work of PGDs86Chapter 6Experimental testing ofYamamoto and Doi theory6.1 Investigation of the J(V) curves under linear sweep ofvoltagesWe start our investigation by subjecting our system to a slow linear sweep of volt-ages ranging from -3 to 3 V at different scanning rates. Below we will show thatunder these conditions the voltammogram obtained via a linear sweep coincide intheir functional form with the predictions of the Yamamoto theory.6.1.1 Voltammograms under linear sweepThe linear sweep of applied voltages produces an J(V) curve shown in Figure 6.1.The figure shows the rectifying mechanism of the PGD. It shows large exponen-tial increase in current for the regime of the forward bias and small currents in theregime of the reverse bias. We note that the regime of forward bias (with exponen-tial increase in currents) doesn’t start until the applied bias is equal to 1.23 V whichis the standard electrode potential of the electrochemical reaction to produce OH– .When the applied voltage is smaller than 1.23 V we observe small negative electriccurrents. The currents increase in magnitude as we increase the negative bias.At∼−2V bias the J(V) curve becomes steeper. The increase in the slope of the87−3 −2 −1 0 1 2 3−0.00100.0010.0020.0030.004V [V] J [A/cm2 ]1.23 V diode breakageFigure 6.1: The resulting IV curves obtained by sweeping through voltages-3 to 3 [V] with the rate of 50 [mV/s].curve is due to the breakage of the diode. In a normal regime the potential dropsat the diode junction prevent the counterions of both sides from mixing. As theapplied electric field becomes larger than the potential drops at the interface, thediode breaks as the counterions are pulled from the opposite side to neutralize theelectrodes. The reduced effective number of counterions results in less H+, OH–being screened. This allows their free diffusion to the interface of the diode andhence, increase in the slope of the J(V) curve is obtained.6.1.2 Regime of forward biasWe also find that experimental evidence of the predictions made by Yamamoto andDoi [77] in the regime of the forward bias. The theory predicts Butler-Volmer likekinetics of water electrolysis at the surfaces of the electrodes. The appeared ionsdiffuse through the system exhibiting exponential currents J/J0 = exp(−∆V )+1.As the forward bias becomes large enough the equation could be further simplified,88Figure 6.2: Lines: the logarithm of the electric current density and the cor-responding data fit for the regime of the forward bias. The findingsconcur with the theoretical prediction of Yamamoto and Doi [77] withJ ∼ exp(−∆V ). Results are obtained by sweeping through voltages -3to 3 [V] with different sweep rates). Circles: the resulting IV curvesobtained by taking terminal values of currents in forward and backwardbias regimes under different magnitudes of applied voltagesince the exponent becomes much larger than 1. Here we find a simple relationshipbetween the current and applied bias J/J0→ exp(−∆V ) when ∆V → ∞.To test the predictions we plotted the log(J) as a function of ∆V (Figure 6.2(lines))which according to Yamamoto theory should be linearly proportional to log(J) ∼∆V . Indeed, we see that starting from 1.23 V the data falls on a linear line rel-atively well. We also find the proportionality constant J0 ≈ e−7 ≈ 1mA/cm2 vialinear regression of the data.89Figure 6.3: Electric current density and the corresponding data fit for theregime of the backward bias. The findings concur with the theoreti-cal prediction of Yamamoto and Doi [77] with J ∼ V 1/2. Results areobtained by sweeping through voltages -3 to 3 [V] with the rate of 50[mV/s]).6.1.3 Regime of backward biasSimilar observations could be made by comparing the predictions of the Yamamotoand Doi [77] theory to experiment in the regime of reverse bias. In this regimethe overpotentials at the electrode surfaces are relatively small, and the diffusionof the H+, OH– are hindered by the floating counterions. As a result the theorypredicts much smaller currents for a given voltage J ∼V 1/2. Indeed by fitting J(V)(Figure 6.3) we find a slope of 1/2 proving the validity of the Yamamoto and Doi[77] in the regime of the backward bias.900 50 100 150 200 250 300time, s2.01.51.00.50.00.51.01.52.0V, [V]0.00030.00020.00010.00000.00010.0002J, [A/cm2]Figure 6.4: The response of a PGD to a step-wise voltage input. Green linesindicate the terminal value of the currents (in the backward and forwardregimes)6.2 PGDs under step voltageHaving understood the behaviour of our systems under a continuous linear sweepwe notice that the curves, for both forward and backward bias remain highly depen-dant on the scanning rates of the sweep. This is unsurprising since application of aever changing potential difference to the electrodes, creates conditions that doesn’tallow the ions to equilibrate. Hence, to analyze the equilibrium currents at theelectrodes J±,as a function of applied bias V± one need to apply a constant voltageand observe evolution of currents as a function of time. Then one could observe aconnection between the currents in the system and the formation of double layer aswell as the ion movement. Based on the assumptions of the Yamamoto theory wewill expect an exponential decay in the observed current, that will equilibrate itselfonto a steady state value of J±. As we expect a rectification of the ionic current weexpect the currents in the forward bias to be much larger than the currents in the91reverse bias, i.e., J+ >> J−.Indeed we observe an exponential decay (increase) in the current value forthe regimes of the forward (backward) biases (Figure 6.4). The curves equilibratewithin τeq ≈ 50−100s suggesting the diffusion constant D≈ L2/τeq ≈ 10−4cm2/sfor the system L ≈ 1 mm wide, which coincides with the theoretical data for theOH– , H+ diffusion 1.After we have performed the experiments with the step biases at different stepbiases ±1 V ,±2 V ,±3 V ,±4 V ,±5 V we aggregate the values J±. We remindthat these values J± represent the terminal values of the current for the regimesof positive and negative bias respectively. We then plot the obtained J currentdensities as a function of applied voltage (Figure 6.2(circles)). As expected, thisgraph, mimics the data for the linear sweep shown in Figure 6.2 with lines as wellas the theoretical predictions of the Yamamoto theory shown in Figure 2.7.Qualitatively Figure 6.2(circles) representing the step bias remains similar tothe regime of a voltage sweep. Quantitatively, however, a difference is observed,between the ionic currents in both forward and backward biases. This differenceoriginates in a slow equilibration of OH– ,H+ from the electrodes to the junctionbetween the gels. When a voltage sweep or step bias are applied, a double layer isformed at the interface between the gel and the electrode. This formation happensrelatively fast within Maxwell-Wagner time of τ = λD/D2 ≈ 1− 100 ns. Then,however, ions require additional 50-100 seconds for slow driven diffusion throughthe sample until they reach the gel junction and recombine to form water. Uponapplication of step bias the system is allowed to equilibrate to equilibrium currentsof J±. When sweep bias is applied (e.g., 50 mV/s) the voltage changes by 50mV/s∗50s = 0.25V within the equilibration time which in its turn changes the observedcurrents. Hence, a careful match against the Yamamoto theory requires a carefulobservation of the equilibrium currents which has been done under the step bias.6.3 Summary• We have verified that basic assumptions of the Yamamoto theory such asequilibration of the current to the steady state coincides with the diffusion1The textbook data of the OH– , H+ was obtained in the regime with no external ions present92times of these of the H+ and OH–• Yamamoto predictions of exponential growth in the forward bias as well as∼V 1/2 dependence in the backward bias are also confirmed• The curves are dependant on the experimental protocol with step bias andslow equilibration is required for careful observations and comparisons againstthe Yamamoto theory93Chapter 7Conclusion and future work7.1 Piezo-ionic effectMolecular dynamics simulations of coarse-grained polyelectrolyte gels have beenperformed to study the Donnan equilibrium and Nernst-Donnan potential at theboundary between gels with different ionizations. We explored several combina-tions of ionizations and different regimes of counterion condensation by varyingthe simulation temperature. The inspiration for the current study comes from thework of Sawahata et al. [64], who argued that the piezoionic (mechanoelectric) ef-fect in gels arises from the stress-induced ionization of its compressed part, thusestablishing a Nernst-Donnan equilibrium. Here, we showed that a reverse method-ology can be used to study this electric potential: two gels with pre-defined ioniza-tion create a lateral pressure gradient that results in a the Nernst-Donnan potentialgenerated at the interface.Depending on the spacing between consecutive backbone macroions (a) andthe Bjerrum length (lB), different types of behavior could be observed. At hightemperatures (small Bjerrum lengths lB < a2 < a1), all counterions are itinerant.They exhibit high osmotic pressure, which creates high elastic pressure of the net-work which expands it and causes high degrees of swelling. Counterions float atrandom positions and the pair energy in the system scales linearly with temper-ature. The high osmotic pressure results in a relatively large break in the localelectroneutrality, and hence large potential gradients occur at the interface between94networks with different ionizations.At low temperatures (high Bjerrum lengths lB > a1,a2) both sides of the gel un-dergo counterion condensation. Positively charged floating counterions on the rightside become bound to negatively charged backbone macroions. The bound coun-terions effectively stop contributing to the osmotic pressure, which causes smallerlevels of swelling of the network on the right hand side. Both parts of the sys-tem experience a decrease in Coulombic energy and an increse in Lennard-Jonesenergy.The Nernst-Donnan potential at the interface scales linearly with temperaturewith the coefficient of proportionality driven by the fraction of concentrations of theuncondensed counterions ln(cfree,2/cfree,1). We also showed that the potential dif-ference can be expressed as a linear function of the lateral pressure. Interestingly,this relationship is robust and hold equally well for weak and strong counterioncondensation. This linear relationship thus provides a molecular level explanationof the piezoionic effect, and calls for verification in experiments.7.2 MD simulation of PE diodeIn order to test the electrostatics of PGDs predicted by a recently introduced PBmodel and to create a modified PB model that more closely reflects the molecu-lar aspects of PGDs, we performed MD simulations of model PGDS at differentelectric field strength, Bjerrum lengths, and degree of ionization. In general, atlow Bjerrum lengths (weak coupling) and low ionization a ”standard” PB modelmatches the MD simulation data rather well and accurately predicts field values atthe boundaries and the peak depth at the interface. With increased strength of elec-trostatic interaction by increasing the Bjerrum length or increasing the ionizationfraction, the standard PB fails to predict the MD profiles due to three additionalphenomena that become relevant: counterion condensation, network collapse atthe interface and field induced gel deformation.At higher Bjerrum lengths lb relative to the inter-backbone ion spacing a (func-tion of ionization degree), we expect counterions to tightly bind to their respec-tive backbone ions and not to contribute to the mobile charge densities. Since thepeaks in the electric field profile occur due to the incomplete charge neutralization,95which is reduced due to the counterion condensation we observe smaller peaks inthe MD simulations. Higher degrees of ionization lead to discrepancies betweenPB and MD results even for low Bjerrum lengths. Due to stronger attraction be-tween the gels in the diode, the collapsed chains at the interface produce a largezone where the polyelectrolyte chains neutralize each other, as well as share theirrespective counterions. This produces a wider interface between the gels with agradual charge density decrease. The ramp in the backbone concentration at theinterface is fundamentally different than the boundary condition used in the PBmodel, where the interface is treated as abrupt. Since in the PB model an externalfield affects only itinerant ions, in the case of MD simulation the polyelectrolytenetwork is allowed to move. This results in the electric field induced deformationof the network, and hence different ion density profiles at the interface and theboundaries of the diode.The counterion condensation coupled with the gel collapse at the interface leadto a much shallower depth of the electric field at the interface. This means thattwo gels sandwiched between electrodes and showing a good rectifying behaviourat low Bjerrum lengths, may show no or much reduced rectifying behaviour ifthe Bjerrum length is increased. It is worth mentioning that although the Bjerrumlength was manipulated artificially in the simulation, it plays an important role inexperimental setups. Depending on the polarization of the solvent, the Bjerrumlength may vary from lB = 7A˙ = 1 σ for a polar solvent like water with εdiel = 78to lB = 35A˙ = 5 σ for less polar solvents like n-propyl alcohol (at T = 40 oC)[2]. Additionally, at higher ionization fractions even small Bjerrum lengths like7A˙ may lead to counterion condensation. Hence changing the solvent, or creatinga PE diode with a high ionization fraction may dramatically change the operatingregime of the diode. As the dip in the electric field profile will decrease, the diodewill operate as a conductor when electric fields of both signs are applied henceceasing to rectify the current.In order to create a model that better mimics the MD results, we modified thePB model to account for counterion condensation as well as the inhomogeneityof the gel backbone density. We found that accounting for a reduction in effectiveionization due to condensation as well as a modification of gel ion density to reflecta ramp up in the backbone ion concentration at the interface produces much shal-96lower electric field profiles. These profiles are in a better agreement with the MDresults.Since the performance of the PGD depends on the size of the potential dipat the interface, accurate prediction of the interface is crucial for proper represen-tation of the diode. The updated PB model informed through MD simulations cannow be combined with a description of the ion currents as described by Yamamotoet al. [77] and should be able to better account for the performance of PGDs indifferent regimes of electrostatic coupling as well as different ionizations.7.3 Experimental investigation of PE diodeIn order to systematically test the predictions of the Yamamoto and Doi [77] theorywe have performed a study of the PGDs under linear sweep and step bias. We havefound that the predictions of the functional forms of I(V) curves of the Yamamototheory hold well. They predict an exponential increase in the regime of the forwardbias as well as the ∼V 1/2 dependency for the regime of the reverse voltage. Bothdependencies were confirmed via a linear sweep of applied biases. Additionally,current equilibration dependency was confirmed via a step voltage bias applied tothe system.7.4 Future workPiezo-ionic sensor:• Future work will entail systematic MD investigation of the piezo-ionic effectin the hydrogel using a fully atomistic model.• For that the PE network will be frozen, and counterions will be allowedto move through the network. Then the affinity of the counterions to thenetwork will be investigated with much better accuracy that coarse-grainedmodel would allow for• This will provide a direct investigation of the diffusion of charged particlesinside of the gel. Furthermore, the MD study will provide a framework uponwhich some of the experimental findings by [61], like diffusion of salt ions97through a charged network could be also explained as a function of affinityof the ions to the network.PE diode:• Yamamoto model will be further stress-tested using accurate 3-electrodesetup measurements to control for the overpotentials upon voltammetric sweepand chronoampermetric measurements.• An electrical impedance spectroscopy (EIS) will be performed to analyzethe structure of the gel and appearance of the double layer at the electrodesurfaces.• Furthermore, depending on the findings provided through experiments Ya-mamoto model could be further enhanced by adding a convection term in theNernst-Plank equations Equation 2.14.• This will summarize the empirical evidence behind the operation of theelectrolysis-based PE diodes and help to optimize their performance.• Additional investigation will be done on non-water electrolysis based diodesbased on ionoelastomer networks ([22, 31]) or electrode-oxidation ([71]) andtheir performance will be compared to the water-electrolysis based ones.98Bibliography[1] Intelligent Hydrogels. Progress in Colloid and Polymer Science, 140:53–62,2013. ISSN 0028-0836. doi:10.1007/978-3-319-01683-2. URLhttp://link.springer.com/10.1007/978-3-319-01683-2. → pages 3, 10[2] G. Akerlof. Dielectric constants of some organic solvent-water mixtures atvarious temperatures. Journal of the American Chemical Society, 54(11):4125–4139, 1932. → pages 28, 96[3] A. Aksimentiev and K. Schulten. Imaging α-hemolysin with moleculardynamics: ionic conductance, osmotic permeability, and the electrostaticpotential map. Biophysical journal, 88(6):3745–3761, 2005. → page 54[4] M. P. Allen et al. Introduction to molecular dynamics simulation.Computational soft matter: from synthetic polymers to proteins, 23:1–28,2004. → pages 42, 44[5] M. Amjadi, K.-U. Kyung, I. Park, and M. Sitti. Stretchable, skin-mountable,and wearable strain sensors and their potential applications: A review.Advanced Functional Materials, 26(11):1678–1698, 2016. ISSN 1616-3028.doi:10.1002/adfm.201504755. URLhttp://dx.doi.org/10.1002/adfm.201504755. → page 1[6] J.-L. Barrat, J.-F. Joanny, and P. Pincus. On the scattering properties ofpolyelectrolyte gels. Journal de Physique II, 2(8):1531–1544, 1992. → page27[7] F. A. Blyakhman, A. P. Safronov, A. Y. Zubarev, T. F. Shklyar, O. A.Dinislamova, and M. T. Lopez-Lopez. Mechanoelectrical transduction in thehydrogel-based biomimetic sensors. Sensors and Actuators A: Physical,248:54–61, sep 2016. → pages 6, 8[8] E. Calo and V. V. Khutoryanskiy. Biomedical applications of hydrogels: Areview of patents and commercial products. European Polymer Journal, 65:99252 – 267, 2015. ISSN 0014-3057.doi:https://doi.org/10.1016/j.eurpolymj.2014.11.024. URLhttp://www.sciencedirect.com/science/article/pii/S0014305714004091. 50Years of European Polymer Journal. → page 1[9] S. Caracciolo, R. G. Edwards, S. J. Ferreira, A. Pelissetto, and A. D. Sokal.Extrapolating monte carlo simulations to infinite volume: Finite-size scalingat ξ/l >> 1. Physical review letters, 74(15):2969, 1995. → page 26[10] O. J. Cayre, S. T. Chang, and O. D. Velev. Polyelectrolyte Diode: NonlinearCurrent Response of a Junction between Aqueous Ionic Gels. Journal of theAmerican Chemical Society, 129(35):10801–10806, 2007. ISSN 0002-7863.doi:10.1021/ja072449z. → pages xiii, 5, 38[11] P.-G. De Gennes. Scaling concepts in polymer physics. Cornell universitypress, 1979. → pages 22, 26[12] Y. Dobashi, M. S. Sarwar, and J. D. W. Madden. MechanoionicTransduction of Solid Polymer Electrolytes and Potential Applications Yuta.MRS Advances, page 1, 2016. doi:10.1557/adv.2016. → pages 4, 5[13] S. Edgecombe, S. Schneider, and P. Linse. Monte Carlo simulations ofdefect-free cross-linked gels in the presence of salt. Macromolecules, 37(26):10089–10100, 2004. ISSN 00249297. doi:10.1021/ma0486391. →page 10[14] A. Erbas and M. Olvera de la Cruz. Energy Conversion in PolyelectrolyteHydrogels. ACS Macro Lett., 4(8):857–861, aug 2015. → pages28, 29, 49, 56, 74[15] A. Erbas and M. Olvera de la Cruz. Energy Conversion in PolyelectrolyteHydrogels. ACS Macro Lett., 4(8):857–861, aug 2015. → pages 10, 11, 12[16] A. Erbas and M. Olvera de la Cruz. Interactions between Polyelectrolyte GelSurfaces. Macromolecules, 49(23):9026–9034, dec 2016. → pages11, 12, 49, 56[17] F. A. Escobedo and J. J. de Pablo. Phase behaviour of model polymericnetworks and gels. Molecular Physics, 90(3):437–444, 1997. ISSN00268976. doi:10.1080/002689797172561. URLhttp://www.tandfonline.com/doi/abs/10.1080/002689797172561.100[18] F. a. Escobedo and J. J. de Pablo. Molecular simulation of polymericnetworks and gels: phase behavior and swelling. Physics Reports, 318:85–112, 1999. ISSN 03701573. doi:10.1016/S0370-1573(99)00012-5.URL http://ac.els-cdn.com/S0370157399000125/1-s2.0-S0370157399000125-main.pdf?{ }tid=ea7c37ee-ba2c-11e4-888a-00000aacb360{&}acdnat=1424566460{ }da9fda98aa12b479d1a60630cd625cc2. → page 10[19] M. Fisher. Comment to ”statistics of long chains with repulsiveinteracations” by j. des cloizeaux. J. Phys. Soc. Jpn., Suppl, 26:44–45, 1969.→ page 26[20] A. Fiumefreddo and M. Utz. Bulk Streaming Potential in Poly(acrylicacid)/Poly(acrylamide) Hydrogels. Macromolecules, 43(13):5814–5819, jun2010. → page 6[21] P. J. Flory. Principles of polymer chemistry. Cornell University Press, 1953.→ pages 22, 26[22] D. Gao and P. S. Lee. Rectifying ionic current with ionoelastomers. Science,367(6479):735–736, 2020. → pages 1, 98[23] F. Gao, F. B. Reitz, and G. H. Pollack. Potentials in anionic polyelectrolytehydrogels. J. Appl. Polym. Sci., 89(5):1319–1321, aug 2003. → pages 6, 8[24] R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L.Seyler, J. Domanski, D. L. Dotson, S. Buchoux, I. M. Kenney, et al.Mdanalysis: a python package for the rapid analysis of molecular dynamicssimulations. Technical report, Los Alamos National Lab.(LANL), LosAlamos, NM (United States), 2019. → page 52[25] H. Guo, T. Kurokawa, M. Takahata, W. Hong, Y. Katsuyama, F. Luo,J. Ahmed, T. Nakajima, T. Nonoyama, and J. P. Gong. QuantitativeObservation of Electric Potential Distribution of Brittle PolyelectrolyteHydrogels Using Microelectrode Technique. J.Chem. Phys., 49(8):3100–3108, apr 2016. → pages 6, 8[26] J. H. Han, K. B. Kim, H. C. Kim, and T. D. Chung. Ionic circuits based onpoly electrolyte diodes on a microchip. Angewandte Chemie - InternationalEdition, 48(21):3830–3833, 2009. ISSN 14337851.doi:10.1002/anie.200900045. → pages 1, 14[27] R. Hockney and J. Eastwood. Computer Simulation Using Particles. Taylorand Francis, jan 1988. doi:10.1201/9781439822050. → page 49101[28] X. Huang, X. Y. Kong, L. Wen, and L. Jiang. Bioinspired Ionic Diodes:From Unipolar to Bipolar. Advanced Functional Materials, 1801079, 2018.ISSN 16163028. doi:10.1002/adfm.201801079. → page 1[29] W. Humphrey, A. Dalke, and K. Schulten. VMD – Visual MolecularDynamics. Journal of Molecular Graphics, 14:33–38, 1996. → page 52[30] E. Jabbari and S. Nozari. Synthesis of acrylic acid hydrogel by g-irradiationcrosslinking of poly (acrylic acid) in aqueous solution. Iran. Polym. J, 8(4):263–270, 1999. → page 21[31] H. J. Kim, B. Chen, Z. Suo, and R. C. Hayward. Ionoelastomer junctionsbetween polymer networks of fixed anions and cations. Science, 367(6479):773–776, 2020. → pages 1, 14, 98[32] D. M. Kochmann, J. S. Amelang, M. I. Espan˜ol, and M. Ortiz. Fromatomistics to the continuum: a mesh-free quasicontinuum formulation basedon local max-ent approximation schemes. PAMM, 11(1):393–394, 2011. →pages xiv, 41[33] P. Kosˇovan, T. Richter, and C. Holm. Modeling of Polyelectrolyte Gels inEquilibrium with Salt Solutions. Macromolecules, 48(20):7698–7708, 2015.ISSN 15205835. doi:10.1021/acs.macromol.5b01428. → pages 3, 10[34] C. Larson, B. Peele, S. Li, S. Robinson, M. Totaro, L. Beccai, B. Mazzolai,and R. Shepherd. Highly stretchable electroluminescent skin for opticalsignaling and tactile sensing. Science, 351(6277):1071–1074, 2016. ISSN10959203. doi:10.1126/science.aac5082. → page 1[35] C. Laschi, B. Mazzolai, and M. Cianchetti. Soft robotics: Technologies andsystems pushing the boundaries of robot abilities. Sci. Robot., 1(1):3690,2016. → page 1[36] E. Lee, J. H. Han, R. Chang, and T. D. Chung. Grand-canonical MonteCarlo simulation study of polyelectrolyte diode. In AIP ConferenceProceedings, volume 1504, pages 1239–1242, 2012. ISBN 9780735411227.doi:10.1063/1.4772152. URLhttp://aip.scitation.org/doi/abs/10.1063/1.4772152. → page 16[37] B. Li, N. Madras, and A. D. Sokal. Critical exponents, hyperscaling, anduniversal amplitude ratios for two-and three-dimensional self-avoidingwalks. Journal of Statistical Physics, 80(3-4):661–754, 1995. → page 26102[38] H. Li, A. Erbas, J. Zwanikken, and M. Olvera de la Cruz. Ionic Conductivityin Polyelectrolyte Hydrogels. Macromolecules, 49(23):9239–9246, dec2016. → pages 16, 17[39] Y. Liu, N. Vrana, P. Cahill, and G. McGuinness. Physically crosslinkedcomposite hydrogels of pva with natural macromolecules: structure,mechanical properties, and endothelial cell compatibility. Journal ofBiomedical Materials Research Part B: Applied Biomaterials: An OfficialJournal of The Society for Biomaterials, The Japanese Society forBiomaterials, and The Australian Society for Biomaterials and the KoreanSociety for Biomaterials, 90(2):492–502, 2009. → page 21[40] A. A. Louis. Beware of density dependent pair potentials. Journal OfPhysics-Condensed Matter, 14(40):9187–9206, 2002. ISSN 0953-8984.doi:10.1088/0953-8984/14/40/311. URLhttp://apps.isiknowledge.com/InboundService.do?product=WOS{&}action=retrieve{&}SrcApp=Papers{&}UT=000179052100012{&}SID=X2hPfKiMI41Oe93ndHk{&}SrcAuth=mekentosj{&}mode=FullRecord{&}customersID=mekentosj{&}DestFail=http://access.isiproducts.com/custom{ }images/wok{ }failed{ }aut. → page 21[41] B. Lovrecek, A. Despic, and J. Bockris. Electrolytic junctions withrectifying properties. The Journal of Physical Chemistry, 63(5):750–751,1959. → page 13[42] S. Mafe´, J. Manzanares, and P. Ramrez. Model for ion transport in bipolarmembranes. Physical Review A, 42(10):6245, 1990. → page 14[43] B. A. Mann. The Swelling Behaviour of Polyelectrolyte Networks. PhDthesis, Johannes Gutenberg-Universita¨t, Mainz, 2005. → pages10, 11, 27, 28, 29, 49[44] B. A. Mann, R. Everaers, C. Holm, and K. Kremer. Scaling inpolyelectrolyte networks. EPL (Europhysics Letters), 67(5):786, 2004. →page 28[45] B. A. Mann, C. Holm, and K. Kremer. Swelling of polyelectrolyte networks.Journal of Chemical Physics, 122(15), 2005. → pages xii, 21, 22, 27, 49[46] B. A. Mann, R. Everaers, C. Holm, and K. Kremer. Scaling inpolyelectrolyte networks. Europhys. Lett., 67(5):786–792, jan 2007. → page49103[47] N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.Mdanalysis: a toolkit for the analysis of molecular dynamics simulations.Journal of computational chemistry, 32(10):2319–2327, 2011. → page 52[48] G. Moy, B. Corry, S. Kuyucak, and S.-H. Chung. Tests of ContinuumTheories as Models of Ion Channels. I. Poisson-Boltzmann Theory versusBrownian Dynamics. Biophysical Journal, 78(5):2349–2363, 5 2000. ISSN0006-3495. doi:10.1016/S0006-3495(00)76780-4. URLhttps://www-sciencedirect-com.ezproxy.library.ubc.ca/science/article/pii/S0006349500767804. → page 15[49] Y. Osada and J.-P. Gong. Soft and wet materials: polymer gels. AdvancedMaterials, 10(11):827–837, 1998. → page 20[50] Y. Osada and A. R. Khokhlov. Polymer gels and networks. 2002. ISBN0203908392. → pages xi, 1, 3, 7, 8, 10[51] S. Plimpton. Fast Parallel Algorithms for Short-Range Molecular Dynamics.Journal of Computational Physics, 117(1):1–19, mar 1995.doi:10.1006/jcph.1995.1039. → page 52[52] S. Plimpton, R. Pollock, and M. Stevens. Particle-mesh ewald and rrespa forparallel molecular dynamics simulations. In PPSC. Citeseer, 1997. → page55[53] K. Prudnikova and M. Utz. Electromechanical Characterization ofPolyelectrolyte Gels by Indentation. Macromolecules, 43(1):511–517, nov2009. → page 9[54] K. Prudnikova and M. Utz. Electromechanical Equilibrium Properties ofPoly(acrylic acid/acrylamide) Hydrogels. J.Chem. Phys., 45(2):1041–1045,jan 2012. → pages 3, 6, 8, 9, 66[55] M. Quesada-Pe´rez, J. Guadalupe Ibarra-Armenta, and A. Martı´n-Molina.Computer simulations of thermo-shrinking polyelectrolyte gels. The Journalof chemical physics, 135(9):094109, 2011. → page 10[56] M. Quesada-Pe´rez, J. Ramos, J. Forcada, and A. Martı´n-Molina. Computersimulations of thermo-sensitive microgels: Quantitative comparison withexperimental swelling data. The Journal of chemical physics, 136(24):244903, 2012.104[57] M. Quesada-Pe´rez, S. Ahualli, and A. Martı´n-Molina. Thermo-responsivegels in the presence of monovalent salt at physiological concentrations: AMonte Carlo simulation study. Journal of Polymer Science, Part B: PolymerPhysics, pages 1403–1411, 2014. ISSN 08876266. doi:10.1002/polb.23576.→ page 10[58] M. Rubinstein and R. Colby. Polymer Physics. OUP Oxford, 2003. ISBN9780198520597. URL https://books.google.ca/books?id=RHksknEQYsYC.→ pages 22, 26, 28[59] M. Rubinstein and R. H. Colby. Polymer physics. 2003. ISBN9780511975691. doi:10.1002/pi.1472. → pages xii, xiii, 3, 23, 24, 26[60] K. Saalwa¨chter, F. Kleinschmidt, and J. U. Sommer. Swellingheterogeneities in end-linked model networks: A combined protonmultiple-quantum NMR and computer simulation study. Macromolecules,37(23):8556–8568, 2004. ISSN 00249297. doi:10.1021/ma048803k. →page 10[61] M. S. Sarwar, Y. Dobashi, E. Glitz, M. Farajollahi, S. Mirabbasi, S. Naficy,G. M. Spinks, and J. D. Madden. Transparent andconformal’piezoionic’touch sensor. In Electroactive Polymer Actuators andDevices (EAPAD) 2015, volume 9430, page 943026. International Societyfor Optics and Photonics, 2015. → pages 4, 5, 6, 97[62] M. S. Sarwar, Y. Dobashi, C. Preston, J. K. Wyss, S. Mirabbasi, and J. D. W.Madden. Bend, stretch, and touch: Locating a finger on an activelydeformed transparent sensor array. Science advances, 3(3):e1602200, 2017.→ pages 1, 2[63] K. Sawahata, M. Hara, H. Yasunaga, and Y. Osada. Electrically controlleddrug delivery system using polyelectrolyte gels. Journal of ControlledRelease, 14(3):253–262, 1990. → page 20[64] K. Sawahata, J. P. Gong, and Y. Osada. Soft and wet touch-sensing systemmade of hydrogel. Macromolecular Rapid Communications, 16(10):713–716, 1995. ISSN 10221336. doi:10.1002/marc.1995.030161002. URLhttp://doi.wiley.com/10.1002/marc.1995.030161002. → pagesxi, 4, 5, 6, 7, 8, 9, 10, 12, 94[65] S. Schneider and P. Linse. Swelling of cross-linked polyelectrolyte gels.European Physical Journal E, 8(5):457–460, 2002. ISSN 12928941.doi:10.1140/epje/i2002-10043-y. → page 10105[66] S. Schneider and P. Linse. Monte Carlo Simulation of Defect-FreeCross-Linked Polyelectrolyte Gels. The Journal of Physical Chemistry B,107:8030–8040, 2003. ISSN 1520-6106. doi:10.1021/jp022336w. URLhttp://pubs.acs.org/doi/abs/10.1021/jp022336w.[67] S. Schneider and P. Linse. Discontinuous volume transitions in cross-linkedpolyelectrolyte gels induced by short-range attractions and strongelectrostatic coupling. Macromolecules, 37(10):3850–3856, 2004. ISSN00249297. doi:10.1021/ma035512n. → page 10[68] A. Sgouros, G. G. Vogiatzis, G. Kritikos, A. Boziki, A. Nikolakopoulou,D. Liveris, and D. Theodorou. Molecular simulations of free and graphitecapped polyethylene films: estimation of the interfacial free energies.Macromolecules, 50(21):8827–8844, 2017. → page 84[69] T. F. Shklyar, A. P. Safronov, O. A. Toropova, G. H. Pollack, and F. A.Blyakhman. Mechanoelectric potentials in synthetic hydrogels: Possiblerelation to cytoskeleton. BIOPHYSICS, 55(6):931–936, mar 2011. → pages6, 8[70] Z. Slouka, M. Prˇibyl, D. Sˇnita, and T. Postler. Transient behavior of anelectrolytic diode. Physical Chemistry Chemical Physics, 9(39):5374–5381,2007. → page 14[71] J.-H. So, H.-J. Koo, M. D. Dickey, and O. D. Velev. Ionic CurrentRectification in Soft-Matter Diodes with Liquid-Metal Electrodes. AdvancedFunctional Materials, 22(3):625–631, 2 2012. ISSN 1616301X.doi:10.1002/adfm.201101967. URLhttp://doi.wiley.com/10.1002/adfm.201101967. → pages 14, 98[72] M. C. Stuart, W. Huck, M. M. Genzer, C. Ober, M. Stamm, G. Sukhorukov,I. Szleifer, V. Tsukruk, U. Marek, F. Winnik, S. Zauscher, I. Luzinov, andS. Minko. Emerging applications of stimuli-responsive polymer materials.Nature Materials, 9(2):101–113, 2010. ISSN 1476-1122.doi:10.1038/nmat2614. URL http://dx.doi.org/10.1038/nmat2614. → page 1[73] A. P. Thompson, S. J. Plimpton, and W. Mattson. General formulation ofpressure and stress tensor for arbitrary many-body interaction potentialsunder periodic boundary conditions. The Journal of chemical physics, 131(15):154107, 2009. → page 54[74] V. Triandafilidi. Coarse-grained molecular dynamics simulation ofpolyelectrolyte gels. https://github.com/bazilevs31/pygels, 2013. → page 52106[75] V. Triandafilidi, S. G. Hatzikiriakos, and J. Rottler. Molecular simulations ofthe piezoionic effect. Soft Matter, 14(30):6222–6229, 2018. ISSN17446848. doi:10.1039/c8sm00939b. → pages 17, 29, 83[76] I. Vlassiouk, S. Smirnov, and Z. Siwy. Nanofluidic Ionic Diodes.Comparison of Analytical and Numerical Solutions. ACS Nano, 2(8):1589–1602, 8 2008. ISSN 1936-0851. doi:10.1021/nn800306u. URLhttp://pubs.acs.org/doi/10.1021/nn800306u. → page 14[77] T. Yamamoto and M. Doi. Electrochemical mechanism of ion currentrectification of polyelectrolyte gel diodes. Nature Communications, 5(May):1–8, 2014. ISSN 20411723. doi:10.1038/ncomms5162. → pagesiv, xiii, xiv, xix, 14, 15, 16, 18, 33, 34, 38, 39, 56, 88, 89, 90, 97[78] Q. Yan and J. J. de Pablo. Monte Carlo simulation of a coarse-grained modelof polyelectrolyte networks. Physical review letters, 91(July):018301, 2003.ISSN 0031-9007. doi:10.1103/PhysRevLett.91.018301. URLhttp://journals.aps.org/prl/pdf/10.1103/PhysRevLett.91.018301. → pages10, 11[79] Q. Yan and J. J. DE PABLO. Monte Carlo Simulation of a Coarse-GrainedModel of Polyelectrolyte Networks. Phys. Rev. Lett., 91(1):018301–4, jul2003. → page 49[80] I.-C. Yeh and M. L. Berkowitz. Ewald summation for systems with slabgeometry. The Journal of Chemical Physics, 111(7):3155–3162, aug 1999.→ page 53[81] J. Zhang, N. Jiang, Z. Dang, T. J. Elder, and A. J. Ragauskas. Oxidation andsulfonation of cellulosics. Cellulose, 15(3):489, 2008. → page 59[82] W. Zhang, X. Zhang, C. Lu, Y. Wang, and Y. Deng. Flexible and transparentpaper-based ionic diode fabricated from oppositely charged microfibrillatedcellulose. Journal of Physical Chemistry C, 116(16):9227–9234, 2012.ISSN 19327447. doi:10.1021/jp301924g. → pages xiii, 14, 38[83] J. Zinn-Justin. Quantum field theory and critical phenomena. Int. Ser.Monogr. Phys., 85:1–996, 1993. → page 26107Appendix ASupporting MaterialsA.1 Validation of the MD modelTo test the usability of the MD model as well as PB code a simplified simulationhas been designed consisting of a simple electrolyte of co-ions and counterions(without the network) placed in a box. The box had periodic boundary conditionsin lateral directions but had reflective sides. Then an electric field Eset has beenapplied to the system and a molecular dynamics simulation has been run. Afterthat a set of several parameters has been compared: Eset , MD profile, analyticalsolution to the linearized PB model as well as a numerical solution for the non-linear PB model both using the in-house code as well as a comercially availablepackage (Comsol) shown in Figure A.1. The figure shows an excellent agreementbetween the curves at least in the limit of weak electric fields. Stronger electricfields, render linearized PB model not applicable and show slight discrepanciesbetween the MD simulations and the solution of the linearized PB model.108Figure A.1: An electrolyte system placed in a box with reflective sides andperiodic boundary conditions in lateral directions and correspondingelectrical field profiles.109