Atomistic simulations of charge transport in photoswitchable organic-graphene hybrids

Photoswitchable self-assembled monolayers (SAMs) in contact with a conductive or semiconductive layer can be used to remotely trigger changes in electrical current using light. In this study, we apply full-atomistic simulations to assess the changes in electronic structure and charge-transport properties of a graphene sheet in contact with an amorphous silica dielectric decorated by an azobenzene SAM. The simulations explicitly account for the structural and electrostatic disorder sourced by the dielectric, which turns out to be weakly affected by photoisomerization and spatially correlated over a length scale of 4–5 nm. Most interestingly, by combining large-scale (tight binding) density functional theory with Kubo–Greenwood quantum transport calculations, we predict that the trans-cis isomerization should induce a shift in surface electrostatic potential by a few tenths of a volt, accompanied by a variation in conductivity by a factor of about 3.


Introduction
Devices in electronics generally feature a single function. While many efforts have been devoted to make them smaller, faster and less energy consumptive, one of the remaining greatest challenge in the field is to impart a multifunctional nature to systems and devices, i.e. to combine electronics with optical, magnetic, or sensing capabilities [1]. Surface coating with self-assembled monolayers (SAMs) constitutes a versatile strategy for introducing in a device additional chemical species that could be exploited for endowing a response to external stimuli. In this context, photochromic molecules are ideal candidates, as light constitutes a fast, non-invasive and low-cost trigger to remotely control the state of the system [2]. Notably, embedding photochromic SAMs in organic or hybrid field-effect transistors opens the way for the fine-tuning of both morphological and electronic properties of the interfaces at stake in the overall response of the device [3][4][5][6]. In the current library of photoswitches, azobenzene derivatives, which exhibit a reversible light-triggered isomerization between a stable trans and a metastable cis state, stand as the most commonly used [7][8][9][10][11][12][13][14][15].
In the following, the focus is on the application of a photoswitchable azobenzene SAM in graphene-based field-effect transistors. Being a 2D atomically thin monolayer, graphene possesses outstanding charge carrier mobility (typically ∼10 3 cm 2 V −1 s −1 , but can be as high as ∼10 5 cm 2 V −1 s −1 [16,17]), while also nowadays available at large scale with well-controlled quality [18,19]. In turn, due to the conformational changes upon irradiation, the two forms of photochromic molecules exert different doping effect, which provides a route to control the carrier type and carrier density of graphene. At the same time, one can expect an improvement of I on /I off ratio, usually rather poor for pristine graphene [20], due to localization effects enabling the development Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. of a transport gap. Besides, unlike the conventional doping via inclusion of the B-and N-atoms or covalent attachment of the functional groups, SAMs should permit higher mobility values since they introduce much smaller perturbations in pi-conjugation within the graphene plane. There has been significant experimental and theoretical work on the design of hybrid materials where electro-active molecules are physisorbed on graphene [21,22], in order to induce surface electrostatic potentials or partial charge transfer at the interface, and in the development of theoretical approaches to predict graphene conductivity [23][24][25]. However, to the best of our knowledge, none of previous studies has gone all the way from an atomistic description of the interfacial microstructure to the modeling of the charge transport properties.
In this work, we analyze the morphological and electronic characteristics of a hybrid assembly of graphene and azobenzene-based SAM (Scheme 1) as a model system for light-responsive field effect transistors. Grafted onto a SiO 2 dielectric, the SAM is in direct contact with the highly sensitive and atomically-thin graphene flake and is expected to perturb its conduction capabilities, thus in turn affecting device performances [26][27][28][29]. The optimization of the interfacial characteristics requires an atomic scale knowledge of the molecular organization in both states of the photoswitch, investigated here using atomistic models. Molecular dynamics (MD) simulations have been widely employed to unravel the structural and dynamical properties of SAMs on SiO 2 [30][31][32][33]. On the other hand, while the detailed description of the photoisomerisation mechanism of azobenzene derivatives calls for high-level quantum-chemical methods [34][35][36][37], such techniques are too expensive for the large surface-SAM systems of interest. Here we resort to a simple classical treatment [38,39] of the photoswitching process to assess the structural organization of SiO 2 /azo-SAM/graphene system in both switching states. The resulting morphologies are fully characterized and used as input for the investigation of charge transport in graphene. To this end, the electrostatic effects induced by the silica-SAM substrate are included in tight-binding (TB) simulations of charge transport within the graphene sheet performed using the Kubo-Greenwood formalism [40].

Methods and computational details
2.1. MD simulations All MD simulations were carried out using the NAMD software [41]. Amorphous silica was described with the Clay force field (FF) [42]. The photochromic SAM molecules were modeled with the General AMBER FF [43] and electrostatic potential fitted (ESP) atomic charges were obtained from quantum-chemical calculations at the ωB97X-D/6-311+G(d, p) level of theory (using the ωB97X-D/6-311 G(d) optimized geometry) using the Gaussian09 software [44]. CASPT2 Ph-N=N-Ph torsional potentials taken from the literature [34] were employed to derive new torsional FF parameters for the ground and excited states S 0 and S 1 of the azobenzene derivatives (figure S1 is available in the supporting information (SI) online at stacks.iop.org/JPMATER/2/ 035001/mmedia) with the methodology described in [45].
From a bulk SiO 2 sample of dimension 124.644×127.328×30 Å 3 , a glass surface was then obtained using an approach similar to that described by Feuston and Garofalini [46,47] following the scheme described in the SI. As a model photoswitchable SAM-forming molecule, we chose the 4-phenyl-azobenzene depicted in Scheme 1, very similar to azobiphenyl thiols studied by Mayor and coworkers on gold surface [8,11,48]. In our case, one of the thiophenyl group has been replaced by a short alkyl chain terminated by a dimethyl-silanolate group to ensure the grafting on silicon oxide [49] and to decouple the photoswitchable azobenzene unit from the surface. Note that, owing to its flexibility and the presence of the two side methyl groups, the alkylsylanol chain is also expected to partially hamper the crystalline packing, which was instead found for more rigid azobiphenyl thiols [11].
Samples of increasing SAM coverage, namely 0.3, 0.8, 1.6, 2.4 and 2.8 molecules nm −2 (composed of 49, 121, 256, 380 and 449 molecules, respectively) were constructed as follows. A photochromic molecule with the azobenzene group in trans conformation was replicated in a regular lattice to generate SAMs of desired densities. Considering the linking of SAM molecules to the surface via a Si-O single bond, the ethyl group of the azoderivatives was removed, and its charge reported on the chromophore oxygen atom. To prepare the surface for the SAM formation, oxygen atoms with coordination number equal to one (located in the first 7 Å of the surface) were removed from the system to provide an adequate number of reactive Si positions. To ensure charge neutrality, the total charge of removed oxygens was redistributed over the atoms of the frozen layer. A specific Lennard-Jones interaction was then defined between silanolate oxygens and silicon atoms from the surface (ε=0.1554 kcal mol −1 , r=1.64 Å), allowing for the two types of atom to be virtually bonded without explicit grafting [30,50,51]. SAM forming molecules were then placed on top of the SiO 2 surface and each sample equilibrated for about 100 ns in the NVT ensemble at 300 K before a production run of 10 ns for characterization of the resulting moprohologies. The thickness is calculated as the projection of the molecular axis, defined as the unit vector joining the silicon atom of the SAM molecule to its terminal carbon, on the surface normal. The tilt angle corresponds to the orientation of the molecular axis with respect to the normal to the surface. The orientational order of the molecules is evaluated from the nematic order parameter P cos , where β is the angle between the molecular axis and the direction of maximum alignment (see [30] for details). Finally, the roughness is calculated as the root mean square deviation of the SAM height (minimum vertical position accessible by a cubic tip with lateral size of 4 Å).
To investigate the effect of a graphene adlayer on the SAM morphology, a graphene conventional 4-atom unit cell of dimensions x=2.456 Å and y=4.256 Å was replicated 51 and 30 times, respectively, to prepare a large layer constituted by 6120 carbon atoms and subsequently placed on top of the equilibrated 2.4 and 2.8 molecules nm −2 SiO 2 /SAM samples. To overcome the mismatch in size of the resulting graphene layer with respect to the dimensions of the SiO 2 surface (and of the simulation box, +0.5% and +0.3% along x and y, respectively), the x and y positions of one carbon atom of the graphene sheet were restrained (harmonic restraint with force constant k=2 kcal mol −1 Å −2 ) thereby allowing the rest of the layer to compress while preventing a possible horizontal shift over the SAM during the simulation. Carbon atoms of graphene were modeled using parameters from Cheng and Steele [52]. Mixed interactions between graphene and SAM molecules were evaluated with Lorentz-Berthelot rules. Note that to preserve the structure of the preformed SAM, the position of the grafting oxygen atom of the molecules was frozen during the graphene deposition and subsequent photoswitching process. After an equilibration stage of 20 ns, a production run of 10 ns was performed.
We adopted an iterative procedure to perform the photoisomerization of the azo-derivatives and obtain the SiO 2 /cis-SAM/graphene system [39,53,54]. Since the purpose of the procedure is to obtain a cis SAM structure, and not to investigate the geometrical changes during the photoisomerization [55], the isomerization process was modeled considering solely the torsional mechanism, using two distinct torsional profiles for the ground (S 0 ) and first excited (S 1 ) states. More specifically, the excitation is modeled by first switching the FF parameterization of the Ph-N=N-Ph torsional potential from the one of S 0 to the one of S 1 . The switch is applied only to molecules in trans conformation, which are thus relaxing in the S 1 potential energy surface for 1 ns (step A, inset of figure S7). The second step (B) consists in turning back the FF parameterization to the ground state potential and let the molecules evolve towards either the cis or trans conformation for 200 ps (details in SI for length of simulations A and B). The procedure is iterated for all unswitched molecules until convergence (whole monolayer switched to cis or 10 consecutive unsuccessful A+B steps). Prior to any characterization of the final cis layer, ESP partial charges of the cis isomer were substituted to the initial set of charges for all the molecules that had switched, while charges of the ones remaining in azo-conformation were left unchanged. With these new atomic charges, the system was equilibrated for 2 ns before a 1 ns production run. The electrostatic potential V(z) was calculated by first computing the charge density as a function of z (the average of the atomic partial charges within slabs parallel to the surface, figure S10), and from it, numerically solving the one-dimensional Poisson equation [55]. To ensure the correctness of this description, MM calculations were checked against DFT to which they correlate satisfactorily (table S1 and figures S4 and S5). The electrostatic potential at each carbon atom, used in the charged transport simulations, was instead computed directly with NAMD by moving a positive test charge at the positions of graphene atoms and taking the mean value as zero of the electrostatic potential scale. Vacuum dielectric permittivity (ε r =0) was used in all electrostatic potential evaluations.

Charge transport simulations
The influence of the SAM-on-silica substrate on the transport properties in the adsorbed graphene monolayers was assessed by performing wave-packet quantum dynamics simulations in the framework of the Kubo-Greenwood formalism [56]. In a drift-diffusion transport regime, the (time and energy dependent) diffusivity, associated to electron wave-packet propagation, saturates to a maximum value, i.e. the diffusion coefficient D max . The (energy dependent) electronic conductivity, σ, relates to this maximum value according to: where ρ(E) is density of states (DOS); the 1/4 factor accounts for 2D transport. As calculations are performed for a periodic system in real space (meaning that only the gamma point is used in the reciprocal space), one has to rely on a graphene sample large enough to achieve the convergence. To tackle this system size, we resort to a model TB Hamiltonian, including up to third nearest-neighbor electronic couplings, parameterized in [25]. The effect of the SAM was mimicked by adding on-site (diagonal) energy disorder to the pristine TB model. The cumulative distribution of energetic disorder was inferred from the in-plane, atomistic, electrostatic potentials computed at the classical level. As discussed below, electrostatic screening effects were either ignored in charge transport simulations, by using the bare electrostatic potentials issued from atomistic simulations (which do not account for intermolecular polarization), or effectively included by assuming a bulk dielectric constant ε r =3. Moreover, we compared two electrostatic disorder realizations characterized by the same magnitude, in terms of standard deviation, but differing in their spatial distribution, namely uncorrelated versus correlated disorder. As noted earlier, the real-space Kubo-Greenwood transport technique requires large system sizes. Therefore, the 51×30 rectangular supercell containing 6120 graphene carbon atoms used in MD simulations was replicated 15 times to obtain a graphene sheet containing 1.377.000 atoms (191.94×188.39 nm 2 ). To avoid spurious spatial correlation effects coming from the replicated images of the electrostatic potential (figure S11), an uncorrelated on-site disorder was introduced for a fraction of randomly chosen sites (10% or 100% of atoms). This disorder was generated by using the inverse sampling technique on the original cumulative distributions of electrostatic potential. Note that the 10% allows for small variations in the shape of the correlation function, and is already sufficient to remove spurious effects arising from the replication of the electrostatic potential of the MD simulation cell, but preserves the value of correlation distance (figure S11). In addition to simulations on the 1.377.000 atoms system, in which the initial wave-packet was completely delocalized, simulations of the dynamic evolution of (initially) localized wave-packet initiated at one carbon atom were performed for a system containing 55.080 carbon atoms (3×3 replication). Temperature effects were accounted for via Fermi-Dirac statistics.
All SAM surfaces exhibit very low roughness, indicating a uniform distribution of molecules on the surface, also at high coverages. Upon increasing the coverage, the thickness of the SAM increases as a result of a gradual change in the layer organization, from a lying-down morphology to a standing-up one (see figure 1 and table 1). This is illustrated by the evolution of the tilt angle θ tilt that decreases from 85°at low coverage, where the molecules lie flat on the surface, to 43°at 2.8 molecules nm −2 (the highest coverage considered in this study). The distribution of cos(θ tilt ) shows that at the intermediate coverage of 1.6 molecules nm −2 the average value of 66°actually results from a broad distribution of the tilt angles (figure 2, light blue full line), and represents an average of tilted and planar groups of molecules (figure 1). More generally, the position of the maximum of the distribution of the tilt angle moves to lower values at increasing coverage, highlighting the transition from planar to vertical alignment.
The orientational order within the SAM was evaluated also through the nematic order parameter 〈P 2 〉, quantifying the overall degree of alignment of the molecular axes (table 1). At very low coverage, 〈P 2 〉 is about 0.3, indicating an almost isotropic distribution of SAM molecules orientations. 〈P 2 〉 remains quite low at intermediate coverage of 1.6 molecules nm −2 , as expected from the disordered organization of the azoderivatives. As the coverage further increases, the molecules straighten up in a more packed and orientationally ordered configuration, resulting in the increase of 〈P 2 〉 to values typical of a nematic liquid crystal phase [57]. 〈P 2 〉 decreases when going from 2.4 to 2.8 molecules nm −2 , suggesting a lower homogeneity of the higher Table 1. Average values of thickness (Å), tilt angle θ tilt (°), nematic order parameter 〈P 2 〉, and rms surface roughness (Å) of the SAM on SiO 2 at different surface coverage (molecules nm −2  density sample, as highlighted also by the snapshot in figure 1 where for the latter coverage, it can be seen the direction of the tilt is not uniform throughout the sample. To get further insight into the relative arrangement of the molecules in the SAM, we calculated the twodimensional radial distribution function (2D RDF) in the x, y plane, which gives a measure of the probability of finding a molecule at a given position with respect to a central reference molecule ( figure S6). The analysis of the 2D RDFs allows estimating the average first neighbor distance, which amounts to about 5 Å and corresponds to the distance between the origin and the first peak of the 2D RDF. However, differently from more standard  SAMs formed by long alkyl chains [30][31][32], no longer-range positional order is present in any of the samples: all SAMs have a liquid-like structure from this point of view.

Effect of a graphene layer on SAM structure
The addition of the graphene layer on top of the SiO 2 /SAM system produces significant changes in the morphology of the 2.4 molecules nm −2 SAM (table 1), namely: (i) a decrease in the thickness of the SAM, (ii) an increase in the value of the tilt angle of its constituting molecules and (iii) a decrease of the roughness. In practice, the strong interaction with graphene makes the SAM to some extent flatter, denser and smoother. Accordingly, the distribution of tilt angles becomes narrower and shifted to higher angles (figure 2, dashed lines). These differences are additionally illustrated by the density profiles across the SAM in absence and presence of graphene (figure 3), which exhibit a sharp peak at the interface with graphene as compared to the smoothly decreasing profile at the interface with vacuum.
The effect of adding graphene on top of the more compact 2.8 molecules nm −2 layer is less significant with respect to the 2.4 molecules nm −2 case, as evidenced by the reduced variations of the tilt angle distribution and the density profiles, the higher density of molecules in this monolayer preventing from larger adjustments of its morphology.

Photoisomerization of the SAM
Starting from the equilibrated all trans samples at 2.4 and 2.8 molecules nm −2 coverage topped with graphene, the photoswitching of the SAMs was performed (see the method section) until convergence at a final isomerization yield of 86% and 82%, respectively (figure S7). Variations in the thickness of azobenzene SAMs of about 4-7 Å when going from trans to cis isomer have been reported by experimental and theoretical studies [58,59], however in our case we do not observe such a marked difference, with a diminution of only about 1 Å (which increases to 2 Å in absence of graphene). This discordance can be attributed to the rather flexible molecular structure of the specific azobenzene derivatives investigated here (Scheme 1), that allows for local conformational deformations, as well as to the relatively disordered packing of the SAMs. To confirm this interpretation, we monitored the change upon isomerization of the orientation for the two molecular fragments at the two ends of the azo group, whose tilt angles with respect to the surface normal are respectively defined by segments 1 and 2 ( figure 4). As shown in figure 4 for a coverage of 2.4 molecules nm −2 , while the segment 1 exhibits a tilt angle on average smaller than the one in the all trans SAM, the upper part of the molecule (segment 2) is more tilted in the cis layer. The combination of these two contrasting effects produces a similar thickness for the all trans and photoisomerised SAMs.
The distribution of tilt angles for segment 2 is rather broad but presents a peak at large angles, indicating that a substantial fraction of SAM molecules have their terminal N-biphenyl moiety parallel and π-stacked to the graphene layer. This feature also translates into a sharper peak in the mass density of the sample at the interface with graphene ( figure 3). Similar but less marked changes are observed at 2.8 molecules nm −2 coverage. We investigated also whether the probability of photoisomerization depends on specific positions on the surface, and in particular on the closeness of already isomerized cis-molecules. Long-range cooperative switching has been experimentally and theoretically demonstrated for azobenzene-based SAMs, and associated with the rigidity of the aromatic backbone, high order, tight packing and π-π interchain interactions [8,9,60], which hinder isomerization unless available space and disorder is created by a neighboring cis-molecule [35]. In our samples, instead, we did not observe any spatial correlation in the bending direction of the switching molecules (figure S8), probably because of the absence of longer-range positional order in the SAM studied here, which is based on rather flexible molecules/anchoring groups.

Variation of the electrostatic potential upon photoisomerisation
Even though the thickness of the SAM is not significantly altered upon photoisomerisation, geometrical modifications occur at the molecular level. Notably, the tilting of the terminal N-biphenyl moiety of the azobenzene molecules is accompanied by a decreased distance between the negatively charged nitrogen atoms and the graphene layer (see figure S9). Such an alteration of the interfacial morphology is expected in turn to bring about a shift in the electrostatic potential exerted by the SAM on the graphene layer. To further investigate this hypothesis, we focused on the 2.4 molecules nm −2 coverage density SAM, which yields the highest structural homogeneity and highest percentage rate of switching. Significant changes indeed appear in the SAM electrostatic potential profile while photoisomerization progresses ( figure 5, left). More specifically, the average electrostatic potential in graphene plane decreases as a linear function of the switch percentage, from 0.47 V for the all trans layer down to −0.06 V in the photoisomerized SAM ( figure 5, right). The origin of this observed negative shift in the electrostatic potential outside the SAM can be traced back to the decrease of the z component of the average dipole moment of the SAM, as a consequence of the random orientation of cis molecular dipole moments in directions mostly parallel to the graphene plane. This in-plane orientation, combined with a slight increase of the molecular tilt angle, leads to the observed decrease of the potential throughout the whole organic layer, upon isomerization. Figure 5 unambiguously reveals a net decrease of the electrostatic potential at the SAM surface upon trans→cis photoisomerization. This negative shift is in line with experimental measurements and periodic calculations for chemically similar SAMs grafted on gold surfaces (shifts of 0.1-0.3 eV [48,61,62]), and consisting in crystalline layers of vertical molecules, while much lower, but always negative shifts were predicted instead for flat lying amorphous layers adsorbed on graphene (0.02-0.04 eV) [21]. Interestingly, our investigation demonstrates that positional order is not strictly necessary for obtaining sizeable voltage shifts, with an overall order of magnitude similar to crystalline SAMs. On the contrary, the orientational order is most likely required since the voltage shifts are determined by the sum of perpendicular components of the molecular dipoles per unit area [63]. Finally, we stress that the predicted average variations in potentials correspond to upper limits, as screening effects are ignored on figure 5. To a first approximation, including a bulk dielectric  constant of 3 to account for screening effects, would yield a change in electrostatic potential of −0.53/3≈ −0.18 V, still a significant number.

Charge transport in the graphene layer
The formation of electron-and hole-rich regions in graphene ('puddles', [64]) is believed to be triggered by the presence of Coulomb impurities associated with the underlying silica surface. The puddles are naturally built-in from our atomistic simulations, as it can be appreciated from the map of the in-plane electrostatic potential shown in figure 6, generated using silica and SAM point charges. Conversely, since we used fixed point charges, electronic polarization and in particular depolarization effects that occur in surfaces where the molecular dipoles are aligned parallel [65,66], are not explicitly accounted for. These effects, on the one hand, produce a screening of the electric charges which is macroscopically represented by a dielectric permittivity ε r higher than unity (close to 3 for azobenzene derivatives [67]) and, on the other, depend nonlinearly on film thickness, with very thin films producing a reduced screening [65].
We then opted for comparing the two limiting cases of ε r =1 (unscreened SAM and SiO 2 charges) and ε r =3, (bulk-like screening) in the charge transport simulations. In the first case we used the electrostatic potential maps for cis and trans as calculated from the MD simulation trajectory, while in the second case all values were scaled down by a factor of 3. In order to characterize quantitatively the spatial inhomogeneity of the electrostatic potential, we plot in figure 7 the distributions of the unscreened electrostatic disorder sampled in the graphene plane, together with the corresponding distance-dependent spatial correlation functions, for the  cis-and trans-azobenzene SAMs: in both cases, we observe a broad Gaussian-like distribution with standard deviation of about 0.58 eV. Taking into account the screening, these values further decrease to 192 mV and 195 mV, respectively, on the order of magnitude of experimental observations of ∼100 mV [68]. As expected from its electrostatic nature, the disorder shows rather long-range spatial correlation with a characteristic scale of about 40 Å for both systems (right panel of figure 7), with (small) deviations in the shape of the curves upon photoswitching. Therefore, photo-isomerization affects both the magnitude of the electrostatic disorder together with its shape, as well as the spatial distribution of the disorder associated with the presence of electronhole puddles. The interplay between these effects is interesting by itself and deserves an additional discussion.
To elucidate the impact of each contribution, we introduced an additional realization of the disorder, namely, spatially uncorrelated (random) disorder, used thereafter for the sake of reference. First, we observe a dramatic impact of the electron-hole puddles on the electronic structure of the graphene layer, as pictured from the DOS in figure 8 (left). More specifically, relative to pristine graphene, the presence of uncorrelated disorder smears out the van Hove singularities but has limited impact in the vicinity of the Dirac point. In contrast, for the same value of standard deviation, the use of correlated disorder, as provided by the FF simulations, entails a much higher DOS in a small energy window around the Dirac point. The presence of the spikes on the DOS Figure 9. Spatial distribution of one-electron wavefunction, 2 y | | (green color, isosurfaces are encompassing 90% of 2 y | | ) at 3.3, 30, and 53 fs, respectively, as extracted from the dynamics of a localized wave-packet for (top) uncorrelated disorder and (bottom) spatially-correlated disorder for the 3×3 replication of the electrostatic potential map of the cis-SAM. The square modulus of the electronic wavefunction 2 y (| | ) was represented using a basis of localized Gaussian functions with standard deviation of 1.42 Å centered on each carbon atom. indicates the localized nature of those states. To provide further insights, we followed the time evolution of an initially localized wave-packet, shown in figure 9 and S12. These simulations revealed that the localization occurs within the puddles as wave-packet is preferentially trapped in the potential well, eventually limiting the charge diffusion at length scale larger than the typical puddle size.
This feature of the electron-hole puddles not only affects the DOS shape but also detrimentally (and negatively) impacts the transport properties, as shown in figure 8 (right). More specifically, we see a reduction of the maximum diffusivity values by an order of magnitude, due to the spatial correlation only. Furthermore, considering screening effects with ε r =3 (center), we observe that the DOS of correlated disorder partly recovers the shape of the pristine graphene, preserving also its spiky character. However, returning to the central problem of the effect of photo-isomerization, the maximum values of diffusion coefficients for the trans and cis conformations are similar no matter which realization of disorder or value of the dielectric constant is considered; therefore, the intrinsic variation in electrostatic disorder upon photo-isomerization results in fairly limited light-induced changes in transport properties.
A more prominent effect is the overall vacuum level shift (VLS) associated with the change in electrostatic potential as the molecules switch from trans to cis. This is shown in figure 10 reporting the conductivity plots (at 300 K) and where the zero of energy is positioned according to the VLS for ε r =1 and ε r =3, respectively. The conductivity versus energy curve mirrors what would be obtained when measuring current versus gate voltage in a transistor setup [69,70]. From figure 10, one can clearly see that over a rather large energy window across the Dirac point, there is significant modulation of the conductivity upon photoswitching. More precisely, the maximum ratio of conductivities (figure 10, right) reaches the value of 2.55 (2.40) at 2.47 eV (0.52 eV) for ε r =1. Furthermore, for ε r =3 we observe the influence of two antagonistic effects, namely, a reduction of the VLS is accompanied by the improvement of transport characteristics. These two effects eventually compensate each other, which yields a similar ratio of conductivities as for ε r =1 with a maximum of 2.65 at 0.31 eV.

Conclusion
In this work, a computational approach encompassing classical MD calculations, DFT and TB electronic structure together with quantum dynamics simulations is applied to model the charge transport properties of a graphene flake in contact with a silica dielectric functionalized with an azobenzene-based photoswitchable SAM. The MD simulations reveal that, irrespective of their conformation, azobenzene molecules self-assemble on silica with limited in-plane positional order and rather high orientational order. The light-induced trans-to-cis isomerization of the grafted azobenzene derivatives is triggered by switching the molecules from their groundto excited-state (torsional) potential energy surface and their relative orientation on the silica substrate tracked as a function of surface coverage. Despite the lack of positional order even at high coverage density, calculations indicate sizeable variations in the average electrostatic potential upon going from the trans to the cis form, as a result of a change in the molecular tilt angle. Besides, the electrostatic potential shows large and correlated spatial fluctuations associated with the silica dielectric, only partly screened by the SAMs. To model large graphene flakes, we resort to a TB model parameterized against DFT electronic structure calculations and map the MD electrostatic potential energy landscape onto corresponding distributions in carbon site energies. Quantum dynamics transport simulations performed in the framework of the Kubo-Greenwood formalism reveal that: (i) the spatially correlated disorder leads to the formation of small puddles trapping charge carriers; and Figure 10. Conductivity in the graphene layer in contact with a cis-and trans-SAM as a function of energy at 300 K, as obtained using correlated disorder for (left) ε r =1 and (center) ε r =3. (Right) ratio in conductivities between trans-and cis-SAMs for ε r =1 (blue line) and ε r =3 (yellow line).
(ii) the trans-to-cis photoconversion is accompanied by a rigid horizontal displacement in the conductivity versus energy curve sparked by an electrostatically driven VLS. In a rather broad energy window around the Fermi energy, we predict a change in electrical conductivity by a factor ∼3 upon photoisomerization. Our results demonstrate that higher sensitivity of the device is hindered by the formation of the electron-hole puddles originating from the defects in SiO 2 . Despite the thickness of ∼11 Å, neither cis-or trans-conformation of the SAM enables an efficient smoothing of the variations in electrostatic potential caused by the substrate. We anticipate that SAM with longer alkyl chains may possess this quality. As an alternative possibility, crystalline substrates may be used to avoid the formation of the puddles. Furthermore, it is also possible to enhance the conductivity variations by combining electrostatic doping coming from the SAM itself with standard electrostatic gating, e.g. via a top gate architecture. Chemical engineering of the SAM forming molecules is also expected to affect theses variations. Finally decreasing the temperature would also increase the sensitivity. We hope these simulations will stimulate experimental work towards the realization of graphene-based lightresponsive devices based on azobenzene SAMs.