Covalent transfer of chemical gradients onto a graphenic surface with 2D and 3D control

Control over the functionalization of graphenic materials is key to enable their full application in electronic and optical technologies. Covalent functionalization strategies have been proposed as an approach to tailor the interfaces’ structure and properties. However, to date, none of the proposed methods allow for a covalent functionalization with control over the grafting density, layer thickness and/or morphology, which are key aspects for fine-tuning the processability and performance of graphenic materials. Here, we show that the no-slip boundary condition at the walls of a continuous flow microfluidic device offers a way to generate controlled chemical gradients onto a graphenic material with 2D and 3D control, a possibility that will allow the sophisticated functionalization of these technologically-relevant materials.


Table of Contents
. In both cases the sealing groove around the microfluidic network is apparent. c, PEEK top layer integrated with the HOPG substrate. d-e, 12 mm x 12 mm HOPG substrate d, before mechanical clamping and e, after the experiment with the mark of the sealing groove created on surface due to the mechanical clamping. Note that the marks generated on the HOPG substrate by the sealing groove allows for mapping out the positions of different zones and areas within these zones, while performing the characterization of the grafted surfaces. Scale bars: 10 mm in a-c, and 2 mm in d-e. Supplementary Fig. 3. Validation of the flow and mass transport in the microfluidic device simulations. a, Velocity profile along x in zone 2 obtained by numerical simulation in the present study (full black line) and from analytical solution 1 (crosses). b, Concentration profile of the NBD diffusion along x in zone 2 obtained by numerical simulation in the present study (full black line), from a 1D numerical solution of the diffusion equation 2 (crosses) and from an analytical solution 3 (circles).
Supplementary Fig. 4. Predictions of aryl radical concentration continuously produced along the channel, obtained by the present numerical model versus an analytical integrated solution. Concentration profile of aryl radical along y (at mid-height, i.e. z = 60 µm, and at mid-width, i.e. x = 0 µm) obtained by numerical simulation in the present study (full black line) and from an analytical solution (crosses) of the integrated reaction rate equation for aryl radical production. In this simulation, only the reaction of aryl radical production was considered, and NBD and KI were assumed to be well-mixed prior to the inlet of the main microfluidic channel and flowed at 100 µL/min (50 µL/min for each of the reactants). In such scenario, the residence times of the fluid elements at the middle of the channel cross-section (z = 60 µm, x = 0 µm) are similar due to the approximately flat velocity profile (Supplementary Fig. 3a), implying similar radical concentrations. Therefore, at that region, diffusion is negligible and the production of radical is the sole mechanism affecting the concentration of aryl radical. This implies that we can compare the concentration profile of aryl radical with the analytical solution of the reaction rate equation, obtained by integration. Supplementary Fig. 5. Numerical simulations. a, Cross-sectional image of xz-plane located in the middle of the main microchannel (i.e. at zone 2) showing the finely meshed domain. b, Concentration map of NBD in the xz-plane at the channel mid-length (i.e., at zone 2). Fig. 6. Velocity map of the flow in the microfluidic channel at mid-height (z = 60 µm) at 600 µL/min. The map shows that the flow becomes fully developed near the beginning of the main microfluidic chamber, meaning that the velocity profiles downstream that region do not change along the length of the chamber. Supplementary Fig. 7. Velocity profile of the developed flow along the width of the main microfluidic channel (x-axis) at mid-height (z = 60 µm) for different flow rates. a-b, Velocity profile in zone 2 along x at mid-height of the channel, obtained from numerical simulations when the flow rate is varied from 10 to 600 µL/min. In panel a, the velocity is expressed in absolute value (mm/s). In panel b, the velocity is expressed in relative terms, by dividing the velocity by its maximum value for each flow rate (relative velocity varies from 0 to 1, where 1 corresponds to the maximum velocity). Fig. 8. Velocity profile of the developed flow along the width of the main microfluidic channel (x-axis) near the HOPG substrate (z = 1 µm) for the different flow rates. a-b, Velocity profile in zone 2 along x near the substrate, obtained from numerical simulations when the flow rate is varied from 10 to 600 µL/min. In panel a, the velocity is expressed in absolute value (mm/s). In panel b, the velocity is expressed in relative terms, by dividing the velocity by its maximum value for each flow rate (relative velocity varies from 0 to 1, where 1 corresponds to the maximum velocity).

Supplementary
Supplementary Fig. 9. Velocity profile of the developed flow along the height of the main microfluidic channel (z-axis) at mid-width (x = 0 µm) for the different flow rates. a-b, Velocity profile in zone 2 along z at mid-width of the channel, obtained from numerical simulations when the flow rate is varied from 10 to 600 µL/min. In panel a, the velocity is expressed in absolute value (mm/s). In panel b, the velocity is expressed in relative terms, by dividing the velocity by its maximum value for each flow rate (relative velocity varies from 0 to 1, where 1 corresponds to the maximum velocity).
Supplementary Fig. 12. Concentration profiles of aryl radical along the main microfluidic chambers' length. a-b, The concentrations of aryl radical, obtained by numerical simulation for a flow rate of 50 µL/min and considering different k1 and different k2, are plotted along the length of the main microchannel (along y) in the middle of its width (x = 0 µm) in close proximity to the HOPG (z = 0 µm). In panel a, the concentration of aryl radical is expressed in absolute value (µM). In panel b, the concentration of aryl radical is expressed in relative terms, by dividing the concentration by its maximum value in each curve (relative concentration varies from 0 to 1, where 1 corresponds to the maximum concentration). The values of k1 and k2 are given in m·s -1 and M -1 ·s -1 , respectively. Supplementary Fig. 13. Concentration profiles of aryl radical along the main microfluidic chambers' length. a-b, The concentrations of aryl radical, obtained by numerical simulation for a flow rate of 50 µL/min and considering k1 = 10 -3 and different k2, are plotted along the length of the main microchannel (along y) in the middle of its width (x = 0 µm) in close proximity to the HOPG (z = 0 µm). In panel a, the concentration of aryl radical is expressed in absolute value (µM). In panel b, the concentration of aryl radical is expressed in relative terms, by dividing the concentration by its maximum value in each curve (relative concentration varies from 0 to 1, where 1 corresponds to the maximum concentration). The values of k1 and k2 are given in m·s -1 and M -1 ·s -1 , respectively. Fig. 14. Rate of aryl radical production, deactivation and grafting along the length of the main microchannel. The reaction rates of aryl radical production (purple blue), deactivation (red) and grafting onto HOPG (cyan), obtained by numerical simulation, are plotted along the length of the main microchannel (along y) in the middle of its width (x = 0 µm) in close proximity to the HOPG (z = 0 µm). The flow rate of each reactant stream was 50 µL/min. The reaction rates shown here were calculated based on equations of aryl radical production, deactivation and grafting, described in the numerical simulation methods, considering k1 = 10 -3 m·s -1 and k2 = 2 x 10 5 M -1 ·s -1 .

Supplementary
Supplementary Fig. 15. Predicted concentration of aryl radicals near the HOPG substrate along the main microfluidic channel's width (x-axis). Three different zones (zone 1-3) along the length of the main microfluidic channel are represented respectively in red, black and blue. The break in the concentration axis is used to clearly show that there are still aryl radicals available at 200-µm far away from the center of the main microfluidic channel (i.e. x = 0). The numerical results show that the aryl radicals can cover the entire width of the main microfluidic channel. Concentration curves calculated based on equations of aryl radical production, deactivation and grafting, described in the numerical simulation methods, considering k1 = 10 -3 m·s -1 and k2 = 2 x 10 5 M -1 ·s -1 . Supplementary Fig. 16. Numerical results of the mass transport in the microfluidic device for a hypothetical case considering free-slip boundary condition at the HOPG substrate. a, Concentration map of the aryl radical concentrations in the entire main microfluidic channel. b, Concentration map of aryl radical at the channel mid-length (i.e. xz-plane, zone 2). c, Concentration map of aryl radical at the channel midlength (i.e. xz-plane, zone 2), considering a very narrow concentration range to show the aryl radical concentration close to the HOPG substrate.  Fig. 5) with their respective line profiles. In panel a, the measured average thickness is 1.1 ± 0.2. However, in this case no clear corral formation is apparent. b, Accordingly, a scratching test is performed and an average thickness of 1.9 ± 0.1 is measured. Note that both analyses result in a different layer thickness since at this flow rate (10 µL/min), the second layer grows on a highly dense monolayer. Scale bars: 500 nm.  Supplementary Table 1. Combinations of reaction rate constants k1 and k2 that were tested to identify those implying aryl radical concentration profiles (in the numerical simulations) that are consistent with the experimentally-observed layer thicknesses. * identifies the combinations that reasonably fit the experimental data, and *** highlights the combination that was the best fit to the experimental data and that was used in the main numerical simulations shown in this work.

Supplementary Note 1: Calculation of diffusion coefficient with DOSY for numerical simulations
The diffusion coefficients were identified by fitting the diffusion profile (the normalized signal intensity as a function of the gradient strength G) at the chemical shift of each signal in the DOSY spectrum with an exponential function using the variant of Stejskal-Tanner equation adapted to the particular pulse sequence used 4 .

Supplementary Note 2: Reaction mechanism
The reaction mechanism affecting the concentration of aryl radicals within the microfluidic channel includes four different main reactions, which involve the formation of the aryl radical and its consumption in the grafting processes, as wells as side reactions leading to the deactivation of the aryl radical 5 . The formation of the aryl radicals occurs when NBD is reduced in a 1-electron process with iodide, resulting in the homolytic breakage of the C-N bond to produce the aryl radical, as well as iodine radicals and N2 gas as by-products. (Supplementary Fig. 27a) The aryl radical formed in the previous step can react with a double bond of the graphitic surface to form the grafted layer. (Supplementary Fig. 27b) Simultaneously, the iodine radicals generated in the first step are highly reactive and can recombine to form molecular iodine. (Supplementary Fig. 27c) The molecular iodine is a well-known scavenger of aryl radicals 6 and reacts rapidly with them to form aryl iodide. (Supplementary Fig. 27d) Since aryl iodide is no longer reactive, the latter step effectively deactivates the radical, as it eliminates its capacity to react with the substrate to form the grafted layer. Therefore, if enough iodine builds up, this reaction will compete with the grafting process and it will be inhibited.

Supplementary Note 3: Calculation of reaction kinetics with UV-Vis spectroscopy for numerical simulations
The method of initial rates was used to determine the rate law of aryl radical formation from NBD and KI to be used in the numerical simulations. The decrease of concentration of NBD was monitored by its absorbance at 260 nm.  Fig. 6). As expected, the velocity of the developed flow along the width of the main microfluidic chamber at mid-height (z = 60 µm) increases significantly with increasing flow rate (Supplementary Fig. 7a).
Moreover, and as expected, the velocity profiles obtained for different flow rates have similar shapes when represented as dimensionless velocity profiles (obtained by dividing each profile by its maximum velocity; Supplementary Fig. 7b). A higher velocity implies a lower residence time during which the reactants can diffuse and react with the HOPG substrate. Importantly, the similar dimensionless velocity profiles for the different flow rates indicate that the velocity of the fluid increases proportionally across the width of the channel, therefore implying proportionally smaller residence times across the width. This indicates that the change in flow rate offers a way to precisely control the residence times, and thus the chemical gradients across the entire width of the microfluidic reactor -key aspects to control the grafting process and the resulting layer formation.
The relation between flow rate and velocity/residence time is similar when the fluid is near the HOPG substrate (z = 1 µm) both in absolute (Supplementary Fig. 8a) and in relative (Supplementary Fig. 8b) terms. However, the absolute velocities were much smaller near the substrate (z = 1 µm) than in the middle of the channel height (z = 60 µm). This is due to the no-slip condition prevailing at the substrate (and all the other surfaces of the microfluidic device) which reduces the velocity of the fluid elements in its vicinity (Supplementary Fig. 9a-b). The no-slip condition is a good approximation of the behavior of most viscous flows near walls, in which the adherence forces between the substrate and the fluid molecules that are closest to it are greater than the cohesion forces between fluid molecules. Because of the no-slip condition, fluid elements in contact with the substrate (and all other surfaces) will have null velocity, and those further away from the substrate/walls will have progressively larger velocities ( Supplementary Figs. 9a-b). Importantly, because the velocity of the flow near the HOPG substrate is much smaller than that in the middle of the channel height (z = 60 µm), the residence times will be much larger, implying larger RD zones. Furthermore, because the height of the device (120 m, z direction) is much smaller than its width (400 m, x direction), the no-slip condition prevailing at the upper and lower surfaces of the microfluidic device leads to velocity profiles along the width (x direction) that are flatter in the central region of the channel (Supplementary Figs. 8a-b). As a consequence, the residence times vary little in that region of the channel, which is particularly convenient to control the reaction time. This shows that the no-slip condition prevailing at the device surfaces is critical to control the substrate functionalization, from the reaction times in the RD zone, to the exact width of substrate to be functionalized.

Supplementary Note 5: Screening of reaction rate constants that fit the experimental data of layer thicknesses
The aryl radical formation, deactivation and grafting processes were modeled based on the equations described in the numerical simulation methods in the main text. Since the rate constants of the grafting reaction (k1) and of the I2 formation (k2) are unknown, we ran simulations for 28 different combinations of k1 and k2 (Supplementary Table 1) to investigate their influence over the aryl radical concentration inside the main microfluidic chamber, and to identify the values fitting best the experimental data of layer thickness.
The value of k1 dictates how fast the grafting reaction occurs. As expected, a large k1 implies fast consumption of aryl radical in the grafting process, leading to a low concentration of aryl radical along the length of the main microfluidic chamber (e.g., blue curve in Supplementary Fig. 10a). Interestingly, when k1 >> 10 0 (for k2 = 10 5 ), the grafting is fast enough to completely consume the aryl radical near the HOPG and, therefore, the concentration of radical becomes limited by the rate of diffusion inside the microfluidic device (which does not depend on the value of k1). This is the case of the simulation with k1 = 10 9 (for k2 = 10 5 , red curve in Supplementary Figs. 10a-b), whose concentration of aryl radical along the length of the main microfluidic chamber is similar (and completely superimposed) to that obtained with k1 = 10 0 (for k2 = 10 5 , blue curve in Supplementary Figs. 10a-b).
On the other hand, when k1 << 10 0 (e.g. 10 -6 , purple curve in Supplementary Fig. 10a), the consumption of aryl radical in the grafting is slow, thus resulting in concentrations of aryl radical that increase along the length of the microfluidic channel (purple curve in Supplementary Figs. 10a-b), despite some deactivation by I2. Moreover, values of k1 of 10 -3 or larger (e.g. green or blue curves in Supplementary Figs. 10a-b), originate concentrations of aryl radical that increase until zone 2, and then decrease until zone 3, as the rate of deactivation starts surpassing the rate of radical formation.
The value of k2 dictates rate of formation of I2, i.e. the species responsible for the deactivation of aryl radical. As expected, a high k2 implies fast deactivation of the aryl radical, and thus a low concentration along the length of the microfluidic chamber (e.g. red curve in Supplementary Figs. 11a-b). More importantly, the value of k2 strongly affects the variation of the radical concentration along the different zones of the microfluidic channel. On one hand, higher values of k2 imply zones 2 and 3 with lower radical concentrations than zone 1, because the fast formation of I2 near the inlet of the main microfluidic chamber induces significant deactivation of aryl radical along the chamber (red curve, Supplementary Fig. 11b). On the other hand, lower values of k2 result in a higher concentration in zone 3 than in zones 1 and 2 because the slow formation of I2 induces little deactivation of aryl radical in the main microfluidic chamber (green curve, Supplementary Fig. 11b). Between these scenarios, there are values of k2 that lead to a higher concentration in zone 2 than in zones 1 and 3, as the formation of I2 occurs throughout the channel (blue curve, Supplementary Fig. 11b).
The k1 and k2 values fitting the experimental data of layer thickness should lead to concentration of aryl radical increasing from zone 1 to zone 2, and decreasing from zone 2 to zone 3, as that would be consistent with layer thicknesses observed in the experiments (Figure 3b). From the 28 combinations tested (Supplementary Table 1), there are multiple combinations that fit this criterion (Supplementary Figs. 10-13). However, most of these combinations only do so marginally, in the sense that they lead to very small increase/decrease between the different zones, that likely would not explain the significant differences in the observed layer thicknesses (see blue and green curves in Supplementary Fig. 12 showing a slight decrease in radical concentration from zone 2 to zone 3). Furthermore, the concentration profiles were found to be quite sensitive to the value of k2, with a change of less than one order of magnitude in k2 resulting in clearly different concentration profiles along the length of the main microfluidic channel (red and blue curves in Supplementary Fig. 13). Our numerical results suggest that the reaction rate constants that best fit the experimental data are k1 = 10 -3 m·s -1 and k2 = 2 x 10 5 M -1 ·s -1 (combination 15 in Supplementary Table   1), as they result in a clear increase in the aryl radical concentration from zone 1 to zone 2, and a similar decrease from zone 2 to zone 3 (blue curve in Supplementary Fig. 13). Moreover, they generate curves of aryl radical concentration along the main reactor (both along its length (y axis, Fig. 2d) and its width (x axis; Fig. 2f)) that are most consistent with the experimentallyobtained thickness of the grafted layer (Fig. 3b). Furthermore, this combination of k1 and k2 was also found to fit the experimental data for other flow rates (Supplementary Fig. 23), which further supports its use to model the dynamics of the HOPG grafting reactions. For this reason, k1 = 10 -3 m·s -1 and k2 = 2 x 10 5 M -1 ·s -1 were used in the numerical simulations reported in the main text of the manuscript. In any case, and despite the obtained consistency between the predicted aryl radical concentrations and the experimentally-observed layer thicknesses, it is important to note that there might be other (k1, k2) pairs that also fit the experimentally-observed layer thicknesses.