# Designing Silicon Brains using LLM: Leveraging ChatGPT for Automated
Description of a Spiking Neuron Array

Mike Tomlinson<EMAIL_ADDRESS>Johns Hopkins University3400 N Charles St,
Baltimore, MD 21218BaltimoreMDUSA21218 ,  Joe Li<EMAIL_ADDRESS>Johns Hopkins
University3400 N Charles St, Baltimore, MD 21218BaltimoreMDUSA21218 and
Andreas Andreou<EMAIL_ADDRESS>Johns Hopkins University3400 N Charles St,
Baltimore, MD 21218BaltimoreMDUSA21218

###### Abstract.

Large language models (LLMs) have made headlines for synthesizing correct-
sounding responses to a variety of prompts, including code generation. In this
paper, we present the prompts used to guide ChatGPT4 to produce a
synthesizable and functional verilog description for the entirety of a
programmable Spiking Neuron Array ASIC. This design flow showcases the current
state of using ChatGPT4 for natural language driven hardware design. The AI-
generated design was verified in simulation using handcrafted testbenches and
has been submitted for fabrication in Skywater 130nm through Tiny Tapeout 5
using an open-source EDA flow.

††conference: ; Jan 25, 2024;

##  1\. Introduction

Over the last three decades, advances in CMOS technology and CAD tools have
led to advances in processor technology that in turn fed research into design
and automation tools that enable the sophisticated System On Chip for general
computing and AI. Verilog and VHDL, both released in the 1980s, have become
standard synthesis tools in digital design. These tools allow the writer to
describe behavioral functionality that can be directly mapped to digital
standard cells and physical layout generation through place and route. Verilog
and VHDL are the schematic entry point into modern CAD tools. Writing and
maintaining code in Verilog and VHDL introduces significant overhead as these
abstract design at at a rather low level. There are a number of projects
trying to address this with a range of adoption and commercial support. These
projects include efforts such as Chisel(Bachrach et al., 2012) and High Level
Synthesis (with specific tools from Cadence, Vivado, and Synopsys). Overall,
the general trend in these methods is to move towards a higher level language
that can then be used to generate VHDL or Verilog.

In November of 2023, OpenAI’s ChatGPT, a LLM captured the attention of users
and business alike because it offered a simple but powerful interface to LLMs
for performing generative AI tasks. This interactive interface to LLMs is
capable of executing a variety of tasks such as writing prose and generating
code. This model has shown to be effective at generating python, albeit with
problems of attention span and adaptability (Kashefi and Mukerji, 2023), (Tian
et al., 2023).

Recent works addressing LLM assisted hardware design include a LLM based
optimization framework that integrates ChatGPT with existing EDA tools. In the
work by (Chang et al., 2023), authors employ LLM tools to implement several
simple modules. For each implementation, power, performance, and area are
compared to modules implemented with ChatGPT alone, Xilinx HLS, and Chisel.
Other research explores a different set of simple blocks, including a simple
microprocessor with a ChatGPT-defined ISA (Blocklove et al., 2023). Similarly,
Yang et al. investigate ChatGPT’s effectiveness for systolic arrays and ML
accelerators (Yang et al., 2023).

###  1.1. Contributions

In this paper we explore the use of generative AI and ChatGPT (version 4 is
used in this work) (OpenAI, 2023) to design a hardware system, namely a
spiking neural network chip, a neuromorphic electronic system (Mead, 1990),
for hardware AI inference (Sanni and Andreou, 2019). Our effort differs from
the current research in the use of LLMs for CAD by focusing on digital spiking
neurons (Cassidy et al., 2013), an unconventional computing architecture, and
by emphasizing complete system design. We document the steps taken to go from
a conversational design description to a functional and synthesizable Verilog
description of a programmable array of spiking neurons. The final AI-generated
HDL has a standardized interface, SPI, and multiple levels of hierarchy. This
work represents one of the first ASICs synthesized entirely from natural
conversational language. By documenting this process, we hope to showcase the
current state of using LLMs as a higher level, conversational, alternative to
handcrafted HDL.

##  2\. Natural Language Hardware Description

This work targets a neuromorphic design with a model trained on a large amount
of public data. Only a small fraction of this data is likely to be related to
neuromorphic engineering and an even smaller fraction to neuromorphic Verilog.
Table 1 gives an idea of the volume of training code available for these
topics. The table lists relevant keywords and the corresponding number of
matching publicly available repositories on Github. There are hundreds of
spiking neuron related repositories, but around 2 orders magnitude fewer
dealing specifically with Verilog.

Table 1. Number of Search Results on Github Search Term | Language | Number of Results
---|---|---
”integrate and fire” | Verilog | 3
”integrate and fire” | Any | 750
”snn” | Verilog | 26
”snn” | Any | 2.4k
”spiking neuron” | Verilog | 6
”spiking neuron” | Any | 478
”spiking” | Verilog | 29
”spiking” | Any | 8.7k
”alu” | Verilog | 1.1k
”accelerator” | Verilog | 339
”cnn” | Verilog | 198
”spi” | Verilog | 530

We start by prompting ChatGPT to generate a leaky integrate and fire (LIF)
neuron. This module is then instantiated in a network module, where 2 layers
of 3 neurons are instantiated in a fully connected fashion. We then prompt for
an SPI peripheral that will be used to program the network parameters.
Finally, we ask ChatGPT to synthesize a top-level module instantiating both
the SPI and network modules. All transcripts and source code can be found on
Github : https://github.com/AndreouLab/tinytapeout_05_chatgpt_snn

###  2.1. The Neuron

The LIF neuron model is one of the more popular spiking neuron models and one
of the simplest. The basic LIF neuron (Cassidy et al., 2013) has one state
variable, the membrane potential, $V_{m}$, and one output, $S$, a boolean
variable indicating if a spike is generated or not. In discrete time, these
variables can be modeled by the following equations:

$V_{m}[t]=\begin{cases}max(I_{in}[t]+V_{m}[t-1]-I_{leak},0),&S[t]=0\\\
0,&S[t]=1\\\ \end{cases}$

$S[t]=V_{m}[t]>V_{th}[t]$

Additional functionality can be added by including a refractory filter. After
a spike occurs, $S[t]$ transitions from 0 to 1; $V_{m}[t]$ will remain at 0
for a fixed period of time. This can be thought of as an additional state
variable.

This discrete time approximation of the LIF neuron model naturally lends
itself to digital circuitry. The membrane potential can be implemented with an
accumulator, the refractory filter, a counter, and the spiking output with a
comparator. With this in mind, we prompted ChatGPT4 with a high level request
for a Verilog module for a spiking LIF neuron:

Prompt: Can you write a verilog module for a spiking leaky integrate and fire
neuron? Please include a multibit current input and a refractory period.
ChatGPT: Sure, I’ll provide a basic Verilog module for a spiking leaky
integrate-and-fire neuron. Here’s a simple version:
(Shown in Table 2.1)

Table 2. Initial and Complete Generated LIF Neuron Verilog

Initial Module |  |  Complete Module
---|---|---
⬇ 1module leaky_integrate_fire_neuron ( 2 input clk, // clock input 3 input reset, // asynchronous reset 4 input [7:0] current, // 8-bit current input 5 output spike // output spike signal 6); 7// Parameters 8parameter THRESHOLD = 8’d255; // spiking threshold 9parameter LEAK_RATE = 8’d1; // decrease by 1 per cycle 10parameter REFRAC_PERIOD = 32; // refractory period in cycles 11 12// Internal states 13reg [7:0] membrane_potential; // current value of the membrane potential 14reg [5:0] refrac_counter = 6’d0; // refractory period counter 15reg in_refrac = 0; // refractory period flag 16 17// On every clock cycle 18always @(posedge clk or posedge reset) begin 19 if (reset) begin 20 membrane_potential <= 8’d0; 21 refrac_counter <= 6’d0; 22 in_refrac <= 0; 23 end else begin 24 if(in_refrac) begin 25 refrac_counter <= refrac_counter - 1’b1; 26 if(refrac_counter == 6’d0) in_refrac <= 0; 27 end else begin 28 // update membrane potential based on current input and leak rate 29 if (membrane_potential + current > THRESHOLD) begin 30 membrane_potential <= 8’d0; 31 spike <= 1; 32 in_refrac <= 1; 33 refrac_counter <= REFRAC_PERIOD; 34 end else if (membrane_potential + current - LEAK_RATE < 8’d0) begin 35 membrane_potential <= 8’d0; 36 end else begin 37 membrane_potential <= membrane_potential + current - LEAK_RATE; 38 end 39 end 40 end 41end 42 43// Reset spike signal after one cycle 44always @(posedge clk or posedge reset) begin 45 if (reset) spike <= 0; 46 else if (spike) spike <= 0; 47end 48endmodule |  |  ⬇ 1module leaky_integrate_fire_neuron ( 2 input clk, // clock input 3 input reset, // asynchronous reset 4 input [7:0] current, // 8-bit current input 5 input [7:0] THRESHOLD, // spiking threshold 6 input [7:0] LEAK_RATE, // decrease by 1 per cycle 7 input [7:0] REFRAC_PERIOD, // refractory period in cycles 8 output reg spike // output spike signal 9); 10// Internal states 11reg [7:0] membrane_potential = 8’d0; 12reg [7:0] refrac_counter = 8’d0; 13reg in_refrac = 0; 14 15// On every clock cycle 16always @(posedge clk or posedge reset) begin 17 if (reset) begin 18 membrane_potential <= 8’d0; 19 refrac_counter <= 8’d0; 20 in_refrac <= 0; 21 spike <= 0; 22 end else begin 23 spike <= 0; // By default, reset the spike signal 24 if(in_refrac) begin 25 refrac_counter <= refrac_counter - 1’b1; 26 if(refrac_counter == 8’d0) in_refrac <= 0; 27 end else begin 28 // Check for potential underflow 29 if (membrane_potential < LEAK_RATE) begin 30 membrane_potential <= current; // Directly set to current 31 end else if (membrane_potential + current < membrane_potential) begin 32 // Potential overflow detected, clamp to THRESHOLD 33 membrane_potential <= THRESHOLD; 34 end else begin 35 membrane_potential <= membrane_potential + current - LEAK_RATE; 36 end 37 38 // Check for spiking 39 if (membrane_potential >= THRESHOLD) begin 40 membrane_potential <= 8’d0; 41 spike <= 1; 42 in_refrac <= 1; 43 refrac_counter <= REFRAC_PERIOD; 44 end 45 end 46 end 47end 48endmodule

After a cursory look at the code, it may seem very impressive! But, after
careful consideration and simulation, a number of problems can be found.

  1. (1)

Line 5: syntax error, spike is declared as a wire (by default), it needs to
explicitly declared as a register

  2. (2)

Lines 27-38: logical error, update logic does not account for overflow

  3. (3)

Lines 34: logical error, faulty logic for checking for an underflow (line 34)

  4. (4)

Lines 18, 45: structural error, spike is multiply driven (used in two
processes)

These problems were all fixable through further prompting. Issues 1 and 4 each
required only one additional prompt to correct. The logical errors, 2 and 3,
were considerably more difficult, requiring several iterations to correctly
handle overflow, underflow, and working with unsigned numbers.

Once the code was functional and synthesizable, we prompted ChatGPT to make
the neuron parameters (threshold, leak rate, and refractory period)
programmable. ChatGPT was able to do this without error, deleting the
parameter declarations and adding input ports to drive the parameter values.
We also increased the bit width of the refractory period from 6 to 8 bits
through prompting. The complete Verilog for this module is shown in Table 2.1.

ChatGPT had little trouble understanding the basic concept of the LIF neuron
and making specific modifications when given a detailed prompt, but struggled
to provide complete modules without syntax errors and showed conceptual
misunderstanding when working with limited precision unsigned numbers. For an
experienced digital designer these problems are relatively easy to fix, but
for a user without the same background, these problems may pose considerable
difficulty. Requiring extensive code modifications heavily limits the ability
for ChatGPT to act as a full natural language abstraction.

###  2.2. Network

For the network, we wanted to implement something relatively simple. We
decided on a 2 layer network with 3 neurons per layer and fully programmable
weights. A diagram of the envisioned network is shown in Fig. 1.

Figure 1. Two layer fully connected neuron network.

Each circle represents one LIF neuron. For a given neuron, the conversion from
spikes to current can be defined as the sum of incoming spikes weighted by the
strength of the synaptic connection. Mathematically, the input current for
each neuron is given by:

$I_{in}[t]=\sum_{i=1}^{N}w_{i}x_{i}[t]$

Where $I_{in}$ is the neuron input current, $x_{i}$ is the incoming spike from
neuron i in the previous layer, $w_{i}$ is the programmed weight for the ith
input to this neuron, and N is the total number of synaptic connections to the
previous layer. To describe this module in Verilog, an experienced designer
might consider using two generate statements, one per layer, and some
combinational logic to calculate the input currents. We opened a new chat,
provided a new copy of the neuron module and then gave the following prompt:

Prompt: Now that we have a working neuron, can you create a new module that
instantiates a network of neurons with 2 layers in a fully connected fashion?
Please connect them with programmable synapses.

The generated code is included in Table LABEL:tab:network. For this module,
there are a number of apparent issues.

  1. (1)

Line 7: syntax error, parameter N is used but never declared

  2. (2)

Lines 7, 9: syntax error, ports, input_currents and weights, are declared as
2d arrays, which is invalid in Verilog (although valid in systemVerilog)

  3. (3)

Line 8: syntax error, similar to the previous issue, spikes is declared as an
unpacked array, which is not allowed for ports in verilog (although valid in
systemVerilog)

  4. (4)

Lines 13-21: syntax error, the instantiation of the neuron is not valid
syntax, no generate variable is used, nor are variables used to index into the
current and spike connections

  5. (5)

Lines 40-48: syntax error, issue 4 also applies to the second set of
instantiations

  6. (6)

Lines 25-37: syntax error, spike variable access is syntactically incorrect.

Issues 1 and 3 were relatively easy to fix, each only requiring one additional
prompt. Issue 2 was more difficult. Eventually, the issue was resolved by
suggesting that ChatGPT flatten the 2D ports into 1D arrays. This required the
prompter to be familiar with Verilog in order to provide the solution,
limiting the use of this process to already experienced engineers.

After flattening the array, ChatGPT created the following code to reassign the
flattened array to an internal 2D variable. This code has a fundamental
problem. An initial block is used when continuous assignment is desired. This
error causes weights to be initialized, but never updated after
initialization, resulting in a module that cannot be programmed. Prompting
ChatGPT to fix this issue caused it to then declare spikes (Line 8) as type
reg, introducing another syntax error (spikes is driven through a port
connection). In fixing one issue, two others were introduced.

⬇

1reg [7:0] weights[2:0][2:0];

2integer i, j;

3initial begin

4 for (i = 0; i < 3; i = i + 1) begin

5 for (j = 0; j < 3; j = j + 1) begin

6 weights[i][j] = weights_flat[(i*3 + j)*8 +: 8];

7 end

8 end

9end

10); \end{lstlisting}

11

12Issue 4 was also difficult to resolve, requiring 7 additional prompts. We
requested ChatGPT explicitly declare each neuron separately and then
concatenate those instantiations back together using a generate statement.
While working on this part of the code, ChatGPT again introduced syntax errors
in other parts of the code. At one point, replacing \textit{begin} in the code
snippet above with an open bracket while keeping the corresponding end
statement, mixing C style coding and Verilog. Similar to the neuron module,
these issues required knowledge of Verilog to identify, debug, and resolve
through prompting, reducing the level of abstraction.

13

14Once the errors were corrected, we used ChatGPT to make design revisions to
the generated code. For this module we added a register file to store the
neuron parameters (shared for all neurons), and the network weights. We first
gave the following prompt:

15

16\hfill \break

17Prompt: \textit{Ok, can you please provide a programmable register file that
stores all the parameters to our network}

18\hfill \break

19

20This prompt produced a reasonable looking register file, that stored the
neuron parameters for each neuron (threshold, leak rate, and refractory
period). The produced code was syntactically correct, but not exactly what we
were asking for. This highlights both an impressive and problematic aspect of
this technology. The previous prompt was ambiguous. It is not clear what
parameters we wanted stored; and rather than asking for clarification, ChatGPT
assumes we want to store independent parameters for each neuron, instead of
storing the weights and a single set of neuron parameters to be shared. This
was an easy clarification to make. The next prompt resulted in an almost
correct module with our desired functionality and only a minor syntax error,
confusing systemVerilog and Verilog coding styles.

21

22\hfill \break

23Prompt: \textit{Can you include weights as a parameter and let’s say both
layers use the same parameters }

24\hfill \break

25

26Again, this prompt was somewhat vague. From the structure of the sentence,
we can infer that the intent is to store independent values for each weight
and a single common copy of the neuron parameters. ChatGPT was able to make
the same inference. Once this module was

27functionally correct, we prompted ChatGPT to merge this functionality into
the network module. This was done without error.

28

29\begin{table*}

30 \caption{Initial Generated Network and SPI modules}

31 \label{tab:network}

32 \begin{tabular}{p{0.45\linewidth} p{0.025\linewidth} p{0.45\linewidth} }

33 \toprule

34 Network Module & & SPI Module \\\

35 \midrule

36 \lstinputlisting[language=verilog]{code/network_v0.v} & &
\lstinputlisting[language=verilog]{code/spi_v0.old.v} \\\

37 \bottomrule

38 \end{tabular}

39\end{table*}

40

41\subsection{Serial Peripheral Interface (SPI)}

42We then prompted ChatGPT to create an SPI peripheral module, allowing us to
program the network parameters and weights over a standard interface.

43

44\hfill \break

45\noindent Prompt: \textit{OK, I want you to now create a SPI interface to
communicate with the network module above.}

46\newline

47

48The generated code is presented in Table \ref{tab:network}, some comments
were removed to reduce the size of the code. This module also has a number of
subtle errors.

49

50\begin{enumerate}

51 \item \textbf{Line 13}: syntax error, use of the SystemVerilog enum feature
which is not supported in Verilog

52 \item \textbf{Lines 27-52}: logical errors in SPI finite state machine

53 \item \textbf{Line 48, 62}: syntax error, write\_enable is multiply driven

54\end{enumerate}

55

56Errors 1 and 3 are similar to errors seen previously. The first is another
mixup with SystemVerilog syntax and the third is a multiply driven net. Both
of these issues were resolved with a single follow-up prompt addressing each.
These consistent mistakes highlight ChatGPT’s lack of familiarity with
Verilog.

57

58Error 2 also showcases a lack of experience with Verilog and with SPI. This
code section confuses one segment and two segment coding styles. Lines 19-25
are typical for a two segment style, but inconsistent with lines 27-52 where
\textit{next\_state} is assigned sequentially, as if \textit{next\_state} is
the state variable in a one segment process. The state logic also has no
dependence on \textit{sclk}, even though it created a port for this signal.
\textit{sclk} is intended to be the SPI clock, which synchronizes the
controller and peripheral. This signal is crucial to writing a correct SPI
state machine.

59

60Additionally, the state machine doesn’t come out of reset. When reset is
applied, state is driven to IDLE and next state is undriven. When reset is
deasserted, the current state will go to whatever \textit{next\_state} was at
the time of applying reset. These issues were all eventually resolved through
subsequent prompts.

61

62% We continued our exploration looking into multiple design alternatives.

63

64% \begin{enumerate}

65% \item list of things we tried

66% \end{enumerate}

67

68% \begin{enumerate}

69% \item Errors encountered on the way

70% \end{enumerate}

71

72

73% Subsequently, after a round of prompting ChatGPT to improve the peripheral
module’s design, it encountered a lapse in consistency. In this instance, it
mistakenly deviated from its own code, removing the reg declarations for the
outputs that it previously established. This glitch resulted in an "illegal
reference to nets" syntax error, shown below:

74% \begin{lstlisting}[language=Verilog]

75% output [3:0] addr_out,

76% output [7:0] data_out,

77% output write_enable

78% \end{lstlisting}

79

80% This again demonstrates that ChatGPT can potentially introduce new
problems while addressing existing ones. We then prompted ChatGPT for more
design enhancements, treating SCLK as a data signal, sampling external signals
using metastability synchronizers, and adding in POCI output. During this
stage, we observed that ChatGPT exhibited minimal errors when tasked with
implementing or modifying specific sections of a pre-existing module. This
held true even when the prompt involved complex tasks, as shown in the
following prompt:

81

82% \hfill \break

83% \noindent

84% Prompt: Make the following changes to the code above: instead of using
sclk like a clk, use it like a normal signal. Add an input clk as system
clock. Furthermore, use a shift register to sample the posedge and neg edge of
sclk.

85% \newline

86

87% ChatGPT demonstrates its ability to accurately handle detailed prompts,
but this can pose as a challenge to less experienced designers. Of course,
even with task-specific prompts, ChatGPT is not completely error-free.
Consider the code below, where ChatGPT attempted to drive the signal
"miso\textunderscore output" in two different always blocks:

88% % \begin{lstlisting}[language=Verilog]

89% % sclk_posedge <= !sclk_prev2 & sclk_prev;

90% % sclk_negedge <= sclk_prev2 & !sclk_prev;

91% % \end{lstlisting}

92

93% % We can see ChatGPT chose logical NOT ’!’ for data assignment when the
best practice is to use bitwise NOT ’ $\sim$’. One may argue this is just a
design choice, but the next problem is a lot more significant, where ChatGPT
attempted to drive the signal "miso\textunderscore output" in two different
always blocks:

94% \begin{lstlisting}[language=Verilog]

95% always @(posedge clk or posedge reset) begin

96% if (reset) begin

97% ...

98% WRITE_EN: begin

99% addr_out <= shift_addr;

100% data_out <= shift_data;

101% write_enable <= 1;

102% miso_output <= {shift_addr, shift_data};

103% state <= IDLE;

104% end

105% ...

106% end

107% end

108% always @(posedge clk) begin

109% if (sclk_negedge) begin

110% miso <= miso_output[11];

111% miso_output <= {miso_output[10:0], 1’b0};

112% end

113% end

114% \end{lstlisting}

115

116% Although this design will pass the compiler and potentially simulation,
it will be rejected by the synthesizer. One plausible explanation for
ChatGPT’s apparent oversight in generating synthesizable code could be
attributed to a lack of labeled samples in its database.

117

118\subsection{Top Module}

119The top module needs to instantiate both the SPI module and the network
module, appropriately connecting internal signals and declaring input and
output ports. In the same chat, the following prompt was used:

120

121\hfill \break

122\noindent Prompt: \textit{Can you create a top file to connect this spi
module with the network module?}

123\newline

124

125The generated code had only one error. The intermediate signals between the
modules were declared as reg instead of wire. This was an easy fix. After one
more prompt, the top module was complete.

##  3\. Tiny Tapeout: Implementation

This design was submitted to TinyTapeout 5, a multi-project die effort for
Skywater 130nm through efabless. The flow is intended to handle most of the
intricacies of digital implementation, exposing only a minimal set of
configuration options in a user level script. This flow uses yosys (Wolf and
Glaser, 2013) for synthesis, which introduced an interesting complication.
Through an iterative debugging process we determined the following code
produced a multiply driven net error.

⬇

1always @(posedge clk or posedge reset) begin

2 if (reset) begin

3 ...

4 for (idx1 = 0; idx1 < 3; idx1 = idx1 + 1) begin

5 FIRST_LAYER_WEIGHTS[idx1] <= 8’d0;

6 ...

7 end

8 end else if (write_enable) begin

9 ...

10 end

11end

12

13always @(*) begin

14 for (idx1 = 0; idx1 < 3; idx1 = idx1 + 1) begin

15 input_currents[idx1] = spikes_in[idx1] ? FIRST_LAYER_WEIGHTS[idx1] : 8’d0;

16 end

17end

18

19Error: : Yosys checks have failed: Encountered check error:

20Warning: Drivers conflicting with a constant 1’0 driver:

21 port Q[2] of cell $procdff$614 ($dff)

After extensive debugging, we found that this error could be resolved by
removing the duplicate loop variables in the code above. Once identified,
simply asking ChatGPT to create new loop variables resolved the issue. Final
simulations using the gate level netlist were conducted through the
TinyTapeout flow. The design occupies 33% of 320um x 200um and is expected to
be fabricated by Summer 2024.

##  4\. Conclusion

This paper explores the use of ChatGPT to convert from natural language to
functionally correct and synthesizable Verilog. We successfully use natural
language entry to generate a complete HDL description of a programmable
spiking neuron array, ready for implementation. It is clear technologies like
ChatGPT have the potential to increase design efficiency, correctly producing
simple modules, quickly generating foundational code from scratch and offering
near instantaneous, accurate modifications of existing code when prompted
detailed instructions. However, the current quality of ChatGPT’s output often
falls short. ChatGPT’s responses tend to frequently include some form of
error, either syntactically or logically. ChatGPT also confidently
demonstrates ignorance of more advanced concepts, leading to potentially
obfuscated bugs, increasing the difficulty of verification. These problems
compound, placing a significant burden on the prompter. If the prompter knows
the solution, they can guide ChatGPT to the answer; but without that
knowledge, it can be difficult to use this technology as a tool for
abstracting Verilog description. Overall, our findings suggest that natural
language to Verilog synthesis has potential; but in its current form, it
leaves much to be desired.

###### Acknowledgements.

This work was supported by NSF Grant 2020624 AccelNet:Accelerating Research on
Neuromorphic Perception, Action, and Cognition through the Telluride Workshop
on Neuromorphic Cognition Engineering, NSF Grant 2332166 RCN-SC: Research
Coordination Network for Neuromorphic Integrated Circuits and NSF Grant
2223725 EFRI BRAID: Using Proto-Object Based Saliency Inspired By Cortical
Local Circuits to Limit the Hypothesis Space for Deep Learning Models.

## References

  * (1)
  * Bachrach et al. (2012) Jonathan Bachrach, Huy Vo, Brian Richards, Yunsup Lee, Andrew Waterman, Rimas Avižienis, John Wawrzynek, and Krste Asanović. 2012.  Chisel: constructing hardware in a Scala embedded language. In _Proceedings of the 49th Annual Design Automation Conference_. ACM, San Francisco California, 1216–1225.  https://doi.org/10.1145/2228360.2228584
  * Blocklove et al. (2023) Jason Blocklove, Siddharth Garg, Ramesh Karri, and Hammond Pearce. 2023.  Chip-Chat: Challenges and Opportunities in Conversational Hardware Design.  http://arxiv.org/abs/2305.13243 arXiv:2305.13243 [cs].
  * Cassidy et al. (2013) Andrew S. Cassidy, Julius Georgiou, and Andreas G. Andreou. 2013.  Design of Silicon Brains in the {nano-}CMOS Era:{} Spiking Neurons, Learning Synapses and Neural Architecture Optimization.  _Neural Networks_ 45 (June 2013), 4–26.  https://doi.org/10.1016/j.neunet.2013.05.011
  * Chang et al. (2023) Kaiyan Chang, Ying Wang, Haimeng Ren, Mengdi Wang, Shengwen Liang, Yinhe Han, Huawei Li, and Xiaowei Li. 2023.  ChipGPT: How far are we from natural language hardware design.  http://arxiv.org/abs/2305.14019 arXiv:2305.14019 [cs].
  * Kashefi and Mukerji (2023) Ali Kashefi and Tapan Mukerji. 2023.  ChatGPT FOR PROGRAMMING NUMERICAL METHODS.  _Journal of Machine Learning for Modeling and Computing_ 4, 2 (2023), 1–74.  https://doi.org/10.1615/JMachLearnModelComput.2023048492
  * Mead (1990) Carver A. Mead. 1990.  Neuromorphic Electronic Systems.  _Proc. IEEE_ 78, 10 (1990), 1629–1636.  https://doi.org/10.1109/5.58356
  * OpenAI (2023) OpenAI. 2023.  GPT-4 Technical Report.
  * Sanni and Andreou (2019) Kayode A. Sanni and Andreas G. Andreou. 2019.  A Historical Perspective on Hardware AI Inference, Charge-Based Computational Circuits and an 8bit Charge-Based Multiply-Add Core in {16nm }FinFET CMOS{}.  _IEEE Journal on Emerging and Selected Topics in Circuits and Systems_ 9, 3 (Sept. 2019), 532–543.  https://doi.org/10.1109/JETCAS.2019.2933795
  * Tian et al. (2023) Haoye Tian, Weiqi Lu, Tsz On Li, Xunzhu Tang, Shing-Chi Cheung, Jacques Klein, and Tegawendé F. Bissyandé. 2023.  Is ChatGPT the Ultimate Programming Assistant – How far is it?  http://arxiv.org/abs/2304.11938 arXiv:2304.11938 [cs].
  * Wolf and Glaser (2013) Clifford Wolf and Johann Glaser. 2013.  Yosys - A Free Verilog Synthesis Suite.  _Proceedings of the 21st Austrian Workshop on Microelectronics(Austrochip)_ (Oct. 2013).
  * Yang et al. (2023) Kaiyuan Yang, Haotian Liu, Yuqin Zhao, and Tiantai Deng. 2023.  A new design approach of hardware implementation through natural language entry.  _IET Collaborative Intelligent Manufacturing_ 5, 4 (2023), e12087.  https://doi.org/10.1049/cim2.12087 _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1049/cim2.12087.


|  CO2-driven diffusiophoresis and water cleaning: similarity solutions for
predicting the exclusion zone in a channel flow†
---|---
|  Suin Shim,a,∗ Mrudhula Baskaran,b Ethan H. Thai,c and Howard A. Stonea,‡
|  We investigate experimentally and theoretically diffusiophoretic separation
of negatively charged particles in a rectangular channel flow, driven by CO2
dissolution from one side-wall. Since the negatively charged particles create
an exclusion zone near the boundary where CO2 is introduced, we model the
problem by applying a shear flow approximation in a two-dimensional
configuration. From the form of the equations we define a similarity variable
to transform the reaction-diffusion equations for CO2 and ions and the
advection-diffusion equation for the particle distribution to ordinary
differential equations. The definition of the similarity variable suggests a
characteristic length scale for the particle exclusion zone. We consider
height-averaged flow behaviors in rectangular channels to rationalize and
connect our experimental observations with the model, by calculating the wall
shear rate as functions of channel dimensions. Our observations and the
theoretical model provide the design parameters such as flow speed, channel
dimensions and CO2 pressure for the in-flow water cleaning systems.

##

††footnotetext: aDepartment of Mechanical and Aerospace Engineering, Princeton
University, Princeton, NJ 08544, USA ††footnotetext: bUnsteady Flow
Diagnostics Laboratory, Institute of Mechanical Engineering, École
Polytechnique Fédérale de Lausanne, 1015, Lausanne, Switzerland††footnotetext:
cDepartment of Electrical Engineering, Princeton University, Princeton, NJ
08544, USA ††footnotetext: † Electronic Supplementary Information (ESI)
available: [details of any supplementary information available should be
included here]. See DOI: 00.0000/00000000.††footnotetext: $\ast$
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
##  1 Introduction

Diffusiophoresis is the spontaneous motion of colloidal particles under a
solute concentration gradient. Since the first recognition 1, 2, the
phenomenon has been widely studied theoretically and experimentally 3, 4, 5,
6, 7, 8, 9, 10, 11, 12, 13, 15, 16, 17, 14, 18, 19. When the solute is an
electrolyte, the net motion of particles has electrophoretic and chemiphoretic
contributions, induced by, respectively, the diffusivity difference(s) among
the ions and the osmotic pressure gradient near the particle surface 5.

Various experimental investigations have been reported for the
diffusiophoresis of micron-sized particles, either in the presence or absence
of liquid flow 12, 13, 22, 14, 18, 19, 20, 21. Recently, it has been
demonstrated that dissolution of gas can drive diffusiophoresis of charged
particles 23, 22, 25, 26, 27, 24. In particular, CO2 was used in several
earlier studies due to its ubiquity and importance in many applications. From
a transport perspective the most important feature of CO2 as a means for
driving diffusiophoresis is that it undergoes rapid dissolution and reaction
in water, creating H+ and HCO${}_{3}^{-}$ ions by the dissociation of H2CO3.
These ions have a large difference in their diffusivities with
$D_{\text{H}^{+}}=9.3\times 10^{-9}$ m2s-1 and
$D_{\text{HCO}_{3}^{-}}=1.2\times 10^{-9}$ m2s-1, leading to a large diffusion
potential and diffusiophoretic mobility of suspended charged particles 23, 22,
25, 26, 27, 24.

As the nature of diffusiophoresis allows moving charged particles with an
imposed directionality, it has been suggested for water cleaning. In our
previous work 26, we demonstrated in a Hele-Shaw geometry that bacterial cells
can be excluded from a region near a CO2 source by CO2-driven
diffusiophoresis. Shin et al.22, Lee et al.21, and Seo et al28, demonstrated
removing particles from water by diffusiophoretic exclusion of particles,
using a CO2 gas chamber or Nafion as the sources of ions, respectively, in a
configuration where the solute gradient is created perpendicular to the liquid
flow. The main mechanism for such a configuration is that either a particle
free zone or an accumulation zone forms near the walls depending on the
diffusiophoretic mobility and the direction of the concentration gradient.
However, the reported experiments use one or two geometrical conditions (in
each study), and we sought to understand what could be achieved in larger
scales with the analogous exclusion zone formation.

Fig. 1:  CO2-driven diffusiophoresis in the channel flow. (a,b) Schematic of
the experimental system. (a) Schematic of the microfluidic channel design. (b)
Zoomed-in view of the channel flow. We define the $x$ coordinate along the
channel starting from the position of the CO2 inlet. (c) Exclusion zone
formation in the flow of negatively charged polystyrene (PS) particles. There
is a finite thickness exclusion zone forming at the lower boundary where CO2
diffuses in. (d) Diffusiophoresis of positively charged amine-modified
polystyrene (a-PS) particles. There is an accumulation of a-PS particles at
the lower (CO2 entering) boundary. Near the top boundary, the particle
concentration decreases, but this is not an exclusion (zero-concentration)
zone. (c,d) The average flow speed is $\langle
u\rangle=100\leavevmode\nobreak\ \mu$m/s. (e) Time evolution of the exclusion
zone in an 1-mm-width channel experiment with the flow speed $\langle
u\rangle=100\leavevmode\nobreak\ \mu$m/s. (f) Visualization of CO2 diffusion
using a pH sensitive dye (Oregon Green 488; Thermo Fisher Scientific). Scale
bars are 500 $\mu$m.

To better understand the transport process in such in-flow water cleaning
systems, and to discuss scale-up possibilities, mathematical models are
needed. A simple increase in the geometry cannot be applied for effective
scale-up of the system due to the nonlinearity of the diffusiophoresis,
possible transitions from laminar to turbulent flows, and other effects (e.g.
gravity) that influence the system as the length scale increases. For example,
CO2 dissolved water is denser than normal water, so a simple increase in the
length scale will cause buoyancy effects that complicate the flow and the
motion of the particles. In such case a different design is necessary.

In this paper, we are interested in a channel flow configuration where the
driving source of diffusiophoretic water cleaning is dissolved CO2, and the
concentration gradient is perpendicular to the flow 22. Through the detailed
investigation of the steady-state flow and the particle exclusion zone using
both theory and experiments, we identify important parameters that directly
affect the scalability of the system and the exclusion zone. In particular, we
apply a shear flow approximation to the near-wall region in the channel flow,
and develop similarity solutions for the distribution of ions and particles.
Then with control experiments and numerical calculations considering the
effect of channel geometry we identify parameters that influences the size of
the exclusion zone, which directly affects the efficiency of the water
cleaning. Our goal is not to provide an optimal solution or the most efficient
channel design. Rather, we try to sort out the important system parameters to
apply to the device design process.

##  2 Formation of the particle exclusion zone in channel flow

We create a concentration gradient of CO2 in a channel flow by placing a CO2
channel and the N2 channel adjacent to the liquid channel. The permeable PDMS
membranes between the channels (thickness $=300\leavevmode\nobreak\ \mu$m) let
the gases dissolve into the liquid channel (see Figure 1(a)). As a result, a
steady state concentration gradient of CO2 is set up in the liquid channel as
described in Figure 1(a,b). In the liquid channel, a particle suspension
(polystyrene; 0.03 vol$\%$ in DI water) flows into the channel at a constant
mean flow velocity. As a result of the CO2 concentration gradient which leads
to the ion (H+ and HCO${}_{3}^{-}$) concentration gradients across the
channel, the charged particles either form an exclusion zone or accumulate
near the CO2-side boundary depending on the surface charge of the particles.
Negatively charged polystyrene (PS) particles form a particle-free zone near
the CO2-side wall as shown in Figure 1(c) while positively charged amine-
modified polystyrene (a-PS) particles accumulate (Figure 1(d)). Details of the
experimental setup is described in the Materials and Methods section. In this
study, we focus on the particle exclusion zone formation near the CO2-side
boundary. Figure 1(e) shows the measured shape of the particle exclusion zone
at different times when negatively charged particles flow into the channel at
a mean speed $\langle u\rangle=100\leavevmode\nobreak\ \mu$m/s. The boundary
of the exclusion zone is plotted versus distance along the channel, which is
measured from the CO2 inlet position (Figure 1(a,b)). We note that some time
after the CO2 stream is turned on, the particle exclusion zone grows then
reach a steady-state as shown in Figure 1(e) (a longer time experiment is
presented in the SI). The nonzero thickness near $x=0$ is due to the CO2
diffusion in the PDMS walls, which makes the influx of CO2 occur from further
upstream. The diffusion of CO2 in the liquid channel can be visualized with a
pH-sensitive dye (Oregon Green 488, Figure 1(f)).

The diffusiophoretic velocity $\mathbf{u}_{p}$ of particles under the ion
concentration gradient $\nabla c_{i}$ is written as

$\mathbf{u}_{p}=\Gamma_{p}\nabla\ln c_{i}=\Gamma_{p}\frac{\nabla
c_{i}}{c_{i}}\leavevmode\nobreak\ ,$ (1)

where $\Gamma_{p}$ is the diffusiophoretic mobility of the particles under a
concentration gradient of a z:z electrolyte 5;

$\Gamma_{p}=\frac{\epsilon}{\mu}\frac{kT}{ze}\left[\beta\zeta-\frac{2kT}{ze}\text{ln}\left(1-\text{tanh}^{2}\frac{ze\zeta}{4kT}\right)\right]\leavevmode\nobreak\
.$ (2)

Fig. 2:  Effect of flow speed and the CO2 pressure on the exclusion zone
formation. (a) Experimental images taken in the 1-mm-width channel at
different flow speeds. For slower particle flow, the exclusion zone thickness
is larger. (b) Measured thickness of exclusion zones obtained at different
flow speeds. (c,d) Effect of CO2 pressure on the exclusion zone formation.
Mean flow velocity is 100 $\mu$m/s. (c) Time evolution of the exclusion zone
thickness at $x=8$ mm for different CO2 pressures. (d) Exclusion zone profile
at $t=300$ s, obtained at different CO2 pressures. Scale bars are 500 $\mu$m.

Here $\epsilon$, $\mu$, $k$, $T$, $e$ and $\zeta$ are, respectively, the
dielectric permittivity of the solution, dynamic viscosity of the solution,
Boltzmann’s constant, absolute temperature and the zeta potential of the
particle. $\displaystyle{\beta=\frac{D_{+}-D_{-}}{D_{+}+D_{-}}}$ is the
diffusivity difference factor, which determines the strength of the local
electric field induced by the difference in diffusivities of the ions.

To obtain initial understanding of the system, we controlled the flow rate in
the liquid channel and varied the CO2 pressure. In the channel flow, the shape
of steady-state exclusion zone can be varied by changing the flow speed
(Figure 2(a,b)). For faster flows, we obtain thinner exclusion zones simply
because the advection and diffusion of CO2 and ions are perpendicular to each
other, and the diffusiophoresis has diffusive feature (note that the
diffusiophoretic mobility has a dimension $\ell^{2}/t$).

In the experiments where CO2 pressure was varied, we observed that there was
not much change in the shape of the exclusion zone. For three different CO2
pressures (2, 10 and 20 psi) applied to the CO2 stream, both time evolution
(Figure 2(c)) and the steady state shape (measured at $t=300$ s, Figure 2(d))
of the exclusion zone were not affected. We interpret that this no-
CO2-pressure-dependence is due to the fast reaction and the logarithmic
feature of diffusiophoresis. From equation (1), we can estimate the particle
velocity by an order of magnitude estimate
$\displaystyle{u_{p}\approx\frac{\Gamma_{p}}{c_{i}}\frac{c_{i}}{w}=\frac{\Gamma_{p}}{w}}$,
where $w$ is the width of the channel over which the ion concentration
gradient is formed. This supports our observation that CO2 pressure did not
affect the shape of particle exclusion zone. The scaling analysis will not
apply to situations where the ion concentration is either very small or very
large, which is considered outside the range of our experiments.

Our initial investigation of the particle exclusion zone with various flow
speeds and CO2 pressures tells us that having slower flow is good for having a
larger exclusion zone, and we only need a moderate amount of CO2 to achieve
such an exclusion zone. However, this is not enough information for scale-up
of the system. For example, slower flow will result in a smaller amount of
clean liquid collected at the outlet, which lowers the system efficiency.
Therefore, we investigate the system in more detail with a mathematical model.

Fig. 3:  Two dimensional model calculation with shear flow approximation. (a)
Schematic of shear flow approximation near the CO2 entering boundary. For
$Pe\gg 1$, the boundary layer approximation can be applied for the diffusion
terms, and we can solve a steady-state, scale invariant problem using the
shear flow approximation. (b) Solutions for equations (8) and (9) for
nondimensional CO2 and ion concentrations. Along the flow – increasing X –,
both CO2 and ions diffuse into the flow in $Y$ direction. Note that the flux
of ions is zero at the CO2-side boundary. (c) Nondimensional particle
distribution along the flow (solutions to equation (12)). Particles are
excluded from the CO2 entering boundary ($Y=0$), and the local accumulation of
particles propagate downstream. (c) Inset: An exclusion zone profile
(arbitrarily defined for the value of $Y$ where $\overline{n}(X,Y)\approx
0.1$) was obtained from the nondimensional particle distribution and plotted
versus $X$. This provides rough estimate of particle free zone profile.

##  3 Steady-state, 2D model with a shear flow approximation

Consider a fully developed flow of a particle suspension in a rectangular
microfluidic channel. If CO2 dissolves into the liquid from one wall as
decribed in Fig. 3(a), we expect to achieve diffusiophoresis of particles
induced by an ion concentration gradient in the channel created by the
reaction 23, 22, 26, 25, 27, 24

$\text{CO${}_{2}$}+\text{H${}_{2}$O}\xrightleftharpoons[k_{r}]{k_{f}}\text{H${}^{+}$}+\text{HCO${}_{3}^{-}$}\leavevmode\nobreak\
.$ (3)

Assume a 2D pressure driven flow with flow velocity u$=u(y)$ ex, then the
reaction-diffusion equation for dissolved CO2 and dissociated ions can be
written as 26, 24

$\displaystyle\frac{\partial c_{c}}{\partial t}+u\frac{\partial
c_{c}}{\partial x}$
$\displaystyle=D_{c}\left(\frac{\partial^{2}c_{c}}{\partial
x^{2}}+\frac{\partial^{2}c_{c}}{\partial
y^{2}}\right)-(k_{f}c_{c}-k_{r}c_{i}^{2}),$ (4) $\displaystyle\frac{\partial
c_{i}}{\partial t}+u\frac{\partial c_{i}}{\partial x}$
$\displaystyle=D_{A}\left(\frac{\partial^{2}c_{i}}{\partial
x^{2}}+\frac{\partial^{2}c_{i}}{\partial
y^{2}}\right)+(k_{f}c_{c}-k_{r}c_{i}^{2}),$ (5)

where $c_{c}$, $c_{i}$ are concentrations of CO2 and ions, and $D_{c}$ and
$D_{A}$ are, respectively, diffusivity of CO2 ($D_{c}=1.9\times 10^{-9}$ m2/s)
and the ambipolar diffusivity of the ions. $k_{f}=0.039$ s-1 and
$k_{r}=9.2\times 10^{4}$ L/mol$\cdot$s are forward and reverse reaction rate
constants 29, 30, 31. The equations can be solved numerically with appropriate
boundary and initial conditions for $u(y)$. For long times we neglect the time
derivative to obtain the steady state concentration field of the ions.

Here we want to seek for solutions for a simplified problem within a small
length scale. Consider CO2 dissolution and dissociation within $x$ small
enough that the concentration distribution of dissolved CO2 does not extend
across the entire channel width. Then the region of interest is a thin layer
near the wall where CO2 is supplied to the liquid phase. Then we can assume
that the liquid flow is simply a shear flow with velocity u$=\dot{\gamma}y$ ex
where $\dot{\gamma}$ is shear rate (Fig. 3(a)-bottom panel), and apply the
boundary layer approximation to further calculations, which is well-known as
the Lévêque approximation in the Chemical Engineering literature. 32.

Now consider steady state reaction-diffusion equations for CO2 and ions in
this near-wall region;

$\displaystyle\dot{\gamma}y\frac{\partial c_{c}}{\partial x}$
$\displaystyle=D_{c}\frac{\partial^{2}c_{c}}{\partial
y^{2}}-(k_{f}c_{c}-k_{r}c_{i}^{2})\leavevmode\nobreak\ ,$ (6)
$\displaystyle\dot{\gamma}y\frac{\partial c_{i}}{\partial x}$
$\displaystyle=D_{A}\frac{\partial^{2}c_{i}}{\partial
y^{2}}+(k_{f}c_{c}-k_{r}c_{i}^{2})\leavevmode\nobreak\ .$ (7)

Boundary conditions are, $c_{c}(x,0)=c_{c}^{\text{sat}}$,
$c_{c}(0,y)=c_{c}(x,\infty)=0$, $\displaystyle{\frac{\partial c_{i}}{\partial
y}\bigg{|}_{(x,0)}=0}$ and $c_{i}(0,y)=c_{i}(x,\infty)=0$, where
$c_{c}^{\text{sat}}$ is the equilibrium concentration of CO2 under the applied
CO2 pressure $p_{c}$ so that $c_{c}^{\text{sat}}=k_{h}p_{c}$ ($k_{h}$ is
Henry’s Law constant). Note that the ions satisfy a zero flux condition at the
CO2-side boundary, and also note that under chemical equilibrium
$k_{f}c_{c}^{\text{sat}}=k_{r}(c_{i}^{\text{sat}})^{2}$, where
$c_{i}^{\text{sat}}$ is the saturation concentration of ions under the applied
gas pressure $p_{c}$. Equations (6) and (7) can be nondimensionalized with
$\overline{c}_{c}=c_{c}/c_{c}^{\text{sat\leavevmode\nobreak\ }}$,
$\overline{c}_{i}=c_{i}/c_{i}^{\text{sat\leavevmode\nobreak\ }}$,
$X=x/\sqrt{D_{c}/\dot{\gamma}}$ and $Y=y/\sqrt{D_{c}/\dot{\gamma}}$. The
nondimensional diffusion-reaction equations are

$\displaystyle Y\frac{\partial\overline{c}_{c}}{\partial X}$
$\displaystyle=\frac{\partial^{2}\overline{c}_{c}}{\partial
Y^{2}}-\frac{k_{f}}{\dot{\gamma}}(\overline{c}_{c}-\overline{c}_{i}^{2})\leavevmode\nobreak\
,$ (8) $\displaystyle Y\frac{\partial\overline{c}_{i}}{\partial X}$
$\displaystyle=\bar{D}_{A}\frac{\partial^{2}\overline{c}_{i}}{\partial
Y^{2}}+\frac{k_{r}c_{i}^{\text{sat}}}{\dot{\gamma}}(\overline{c}_{c}-\overline{c}_{i}^{2})\leavevmode\nobreak\
,$ (9)

where $\bar{D}_{A}=D_{A}/D_{c}$ is the nondimensional ambipolar diffusivity.
The nondimensional boundary conditions are then, $\overline{c}_{c}(X,0)=1$,
$\overline{c}_{c}(0,Y)=\overline{c}_{c}(X,\infty)=0$,
$\displaystyle{\frac{\partial\overline{c}_{i}}{\partial Y}\bigg{|}_{(X,0)}=0}$
and $\overline{c}_{i}(0,Y)=\overline{c}_{i}(X,\infty)=0$. We note that
$k_{f}/\dot{\gamma}$ is dimensionless, and for $\dot{\gamma}=2$ s-1,
$k_{f}/\dot{\gamma}\approx 0.02\ll 1$. Equations (8) and (9) were solved with
Matlab and plotted versus $Y$ for different $X$ values in Figure 3(b). The
plots show that both CO2 and ions diffuse into the channel (toward larger $Y$)
as the liquid flows downstream (increasing $X$). Note that the CO2 has a
constant pressure (and so the concentration) condition at the boundary and the
solutions for ions show zero-flux condition at the boundary. Using the
concentration of ions and the diffusiophoretic velocity of particles, now we
can solve for the particle distribution in the channel also under the shear
flow approximation.

The 2D unsteady advection-diffusion equation for particles of concentration
$n$ is

$\frac{\partial n}{\partial t}+u\frac{\partial n}{\partial
x}+\frac{\partial(u_{p,x}n)}{\partial x}+\frac{\partial(u_{p,y}n)}{\partial
y}=D_{p}\left(\frac{\partial^{2}n}{\partial
x^{2}}+\frac{\partial^{2}n}{\partial y^{2}}\right)\leavevmode\nobreak\ ,$ (10)

where $u_{p,x}$ and $u_{p,y}$ are, respectively, the $x$ and $y$ components of
diffusiophoretic velocity $u_{p}$. For small $x$ and $y$, we can apply the
shear flow approximation to the particle equation. Note that under the
boundary layer approximation, we can neglect the diffusiophoresis in the $x$
direction for
$\displaystyle{\frac{\Gamma_{p}}{x}<\dot{\gamma}y\leavevmode\nobreak\
\Rightarrow\leavevmode\nobreak\
\frac{\Gamma_{p}}{\dot{\gamma}(D_{c}/{\dot{\gamma}})}=\frac{\Gamma_{p}}{D_{c}}<1}$.
This is valid for most experimental situations so we keep the condition for
further calculations. Then we obtain a steady-state equation for particles,

$\dot{\gamma}y\frac{\partial n}{\partial x}+\Gamma_{p}\frac{\partial}{\partial
y}\left(n\frac{\partial\ln c_{i}}{\partial
y}\right)=D_{p}\frac{\partial^{2}n}{\partial y^{2}}\leavevmode\nobreak\ .$
(11)

Boundary conditions are, $n(0,y)=n(X,\infty)=n_{0}$ and
$\displaystyle{-\Gamma_{p}\frac{\partial\ln c_{i}}{\partial
y}n+D_{p}\frac{\partial n}{\partial y}=0}$ at $y=0$, where $n_{0}$ is the
unperturbed particle concentration.

The equation can be nondimensionalized with $\overline{n}=n/n_{0}$,
$X=x/\sqrt{D_{c}/\dot{\gamma}}$ and $Y=y/\sqrt{D_{c}/\dot{\gamma}}$, and thus

$Y\frac{\partial\overline{n}}{\partial
X}+\bar{\Gamma}_{p}\frac{\partial}{\partial
Y}\left(\overline{n}\leavevmode\nobreak\
\frac{\partial\ln\overline{c}_{i}}{\partial
Y}\right)=\bar{D}_{p}\frac{\partial^{2}\overline{n}}{\partial
Y^{2}}\leavevmode\nobreak\ ,$ (12)

where $\bar{\Gamma}_{p}=\Gamma_{p}/D_{c}$ and $\bar{D}_{p}=D_{p}/D_{c}$.
Equation (12) is solved and the solution is plotted versus $Y$ for different
positions along the channel in Figure 3(c). Note that the solutions show zero
particle concentration region near $Y=0$ which indicates the particle
exclusion zone. We calculated the value of $Y$ where $\overline{n}(X,Y)=0.1$,
then plotted this location versus $X$ (inset of Figure 3(c)). This plot shows
an approximate shape of the steady-state particle exclusion zone obtained from
the model calculation. For the shear rate $\dot{\gamma}=5$ s-1, the
nondimensional thickness of exclusion zone $Y_{0.1}\approx 4$ yields a
dimensional length scale is $y_{0.1}\approx 80\leavevmode\nobreak\ \mu$m,
which provides a reasonable estimate of a typical length scale of exclusion
zone thickness.

Performing numerical calculation of the partial differential equations (PDEs)
provides a good estimate of the particle exclusion zone, but we can further
simplify our equations. From the form of the equation, we can deduce a
similarity variable and analyze the system with a set of ordinary differential
equations (ODEs).

##  4 Similarity solutions for predicting particle exclusion zone

Fig. 4:  Similarity solutions for concentrations of the CO2, ions and the
particles. (a) Similarity solution for CO2 concentration plotted versus
$\eta$. Inset: CO2 concentration plotted versus $Y$ for different values of
$X$. (b) Ion concentration $f_{i}(\eta)$ plotted versus $\eta$. Inset: ion
concentration $\overline{c}_{i}$ is plotted versus $Y$ for different $X$’s.
Note that the inset graphs in (a) and (b) are solutions to approximated
equations, and so they show different solutions from Figure 3. (c, d)
Solutions of equation (eq) are plotted for different values of
$\bar{\Gamma}_{p}$. (c) Positively charged particles ($\bar{\Gamma}_{p}>0$)
tend to accumulate more as the mobility increases near $\eta\rightarrow 0$.
(d) Negatively charged particles ($\bar{\Gamma}_{p}>0$) move away from the
boundary ($\eta=0$), and the larger the mobility is, more particles are
cleared away from the CO2 boundary. There is a finite value of $\eta$ for each
negative $\bar{\Gamma}_{p}$ which indicate the thickness of particle free zone
in terms of similarity variable.

From the structure of equation (8), we can consider a similarity variable
$\displaystyle{\eta=\frac{Y}{g(X)}}$ and find a similarity solution
$f_{c}(\eta)=\overline{c}_{c}(X,Y)$. Substituting $\overline{c}_{c}$ with
$f_{c}(\eta)$ in equation (8), we get

$-\eta^{2}\frac{df_{c}}{d\eta}g^{\prime}g^{2}=\frac{d^{2}f_{c}}{d\eta^{2}}-\frac{k_{f}g^{2}}{\dot{\gamma}}(f_{c}-\overline{c}_{i}^{2})\leavevmode\nobreak\
.$ (13)

Finding $g(X)$ that satisfies $g^{\prime}g^{2}=1$ gives
$g(X)=(3X)^{\frac{1}{3}}$, and thus equation (13) becomes

$-\eta^{2}\frac{df_{c}}{d\eta}=\frac{d^{2}f_{c}}{d\eta^{2}}-\frac{k_{f}}{\dot{\gamma}}(f_{c}-\overline{c}_{i}^{2})(3X)^{\frac{2}{3}}\leavevmode\nobreak\
.$ (14)

Since $\displaystyle{\frac{k_{f}}{\dot{\gamma}}X^{\frac{2}{3}}\ll 1}$ for
small $X$ (note that $k_{f}/\dot{\gamma}\approx 0.02$ for $\dot{\gamma}=2$
s-1), the reaction term can be neglected. Therefore, we solve an ODE

$-\eta^{2}\frac{df_{c}}{d\eta}=\frac{d^{2}f_{c}}{d\eta^{2}}\leavevmode\nobreak\
,$ (15)

with the boundary conditions $f_{c}(0)=1$ and $f_{c}(\infty)=0$. This is a
similarity transform of the well-known Lévêque problem, and the solution to
equation (15) is

$f_{c}(\eta)=\frac{\int_{\eta}^{\infty}e^{-s^{3}/3}ds}{\int_{0}^{\infty}e^{-s^{3}/3}ds}=\frac{\Gamma\left(\frac{1}{3},\frac{\eta^{3}}{3}\right)}{\Gamma\left(\frac{1}{3}\right)}\leavevmode\nobreak\
.$ (16)

Here $\Gamma(\cdot)$ is the Gamma function and $\Gamma(\cdot,\cdot)$ is the
incomplete Gamma function, i.e.
$\displaystyle{\Gamma\left(\frac{1}{3},\frac{\eta^{3}}{3}\right)=\int_{\frac{\eta^{3}}{3}}^{\infty}t^{-\frac{2}{3}}e^{-t}dt}$.
The solution (16) is plotted in Figure 4(a). The inset of Figure 4(a) shows
the solutions to the Lévêque problem plotted versus $Y$ at different values of
$X$. The plots show the diffusion of CO2 into the channel starting from the
boundary ($\eta=Y=0$).

The same analysis can be applied to the concentration of ions near the
boundary. For equation (9), we can find a similarity solution $f_{i}(\eta)$ by
defining
$\displaystyle{\overline{c}_{i}(X,Y)=\frac{k_{r}c_{i}^{\text{sat}}}{\dot{\gamma}}(3X)^{\frac{2}{3}}f_{i}(\eta)}$
with $\displaystyle{\eta=\frac{Y}{(3X)^{\frac{1}{3}}}}$. Substituting this
into equation (9), we obtain

$2\eta
f_{i}-\eta^{2}\frac{df_{i}}{d\eta}=\bar{D}_{A}\frac{d^{2}f_{i}}{d\eta^{2}}+\left(f_{c}-\left(\frac{k_{r}c_{i}^{\text{sat}}}{\dot{\gamma}}\right)^{2}(3X)^{\frac{4}{3}}f_{i}^{2}\right)\leavevmode\nobreak\
,$ (17)

where the term
$\displaystyle{\left(\frac{k_{r}c_{i}^{\text{sat}}}{\dot{\gamma}}\right)^{2}(3X)^{\frac{4}{3}}\ll
1}$ for small $X$. Therefore,

$2\eta
f_{i}-\eta^{2}\frac{df_{i}}{d\eta}=\bar{D}_{A}\frac{d^{2}f_{i}}{d\eta^{2}}+f_{c}\leavevmode\nobreak\
,$ (18)

and the boundary conditions are $f_{i}^{\prime}(0)=0$ and $f_{i}(\infty)=0$.
The ion concentration satisfies no-flux condition at the CO2-side boundary.
This equation shows that the ion concentration depends on the CO2
concentration by the forcing term $f_{c}$. Equation (18) is solved and plotted
in Figure 4(b). The inset of Figure 4(b) is the plot of
$\displaystyle{\overline{c}_{i}(X,Y)=\frac{k_{r}c_{i}^{\text{sat}}}{\dot{\gamma}}(3X)^{\frac{2}{3}}f_{i}(\eta)}$
versus $Y$ at different values of $X$. For increasing $X$, the ions created by
the reaction diffuse into the channel.

Finally, we can consider the particle distribution near the CO2 boundary where
CO2 enters. Define $N(\eta)=\overline{n}(X,Y)$ and substitute into equation
(12) we obtain

$-\eta^{2}\frac{dN}{d\eta}+\bar{\Gamma}_{p}\frac{dN}{d\eta}\frac{d\ln
f_{i}}{d\eta}=\bar{D}_{p}\frac{d^{2}N}{d\eta^{2}}-\bar{\Gamma}_{p}N\frac{d^{2}\ln
f_{i}}{d\eta^{2}}\leavevmode\nobreak\ .$ (19)

The boundary conditions are, $\displaystyle{-\bar{\Gamma}_{p}N\frac{d\ln
f_{i}}{d\eta}+\bar{D}_{p}\frac{dN}{d\eta}=0}$ at $\eta=0$ and $N(\infty)=1$.
Under typical experimental conditions, $\bar{D}_{p}\ll\bar{\Gamma}_{p}$, and
thus we can neglect the particle diffusion. As a result, we solve a first-
order ODE for $N(\eta)$,

$-\eta^{2}\frac{dN}{d\eta}+\bar{\Gamma}_{p}\frac{d}{d\eta}\left(N\frac{d\ln
f_{i}}{d\eta}\right)=0\leavevmode\nobreak\ ,$ (20)

with a boundary condition $N(\infty)=1$. Both for positively charged
($\bar{\Gamma}_{p}>0$) and negatively charged particles
($\bar{\Gamma}_{p}<0$), the equation (20) is solved and plotted versus $\eta$
in Figure 4(c,d) with five different magnitudes of the mobility
$\bar{\Gamma}p$. Solutions to the full equation (19) are plotted in the SI.
For positively charged particles (Figure 4(c)), there is accumulation
(increase in particle concentration) near the boundary $\eta=0$, and the
particle accumulation is enhanced for larger mobilities. For negatively
charged particles (Figure 4(d)), we note that for each value of
diffusiophoretic mobility, there is a critical value of $\eta$ below which the
particle concentrations are effectively zero. We call this cutoff value
$\eta_{c}$. In our experiments, the diffusiophoretic mobility of the
negatively charged polystyrene particles corresponds to
$\bar{\Gamma}_{p}=-0.4$, and the critical value $\eta_{c}=0.78$.

Fig. 5:  Comparison between experiments and equation (22) for different flow
speeds. Experiments are done in a 1-mm-width and a 100-micron-height channel.
With the choice of $\alpha=0.45$, equation (22) fits the experimental
measurements under different conditions. Fig. 6:  Calculation of flow
velocities and shear rates for different channel dimensions. (a) Schematic of
the rectangular channel. (b) Height-averaged velocity profiles
($\bar{u}(y)/\langle u\rangle$) plotted versus $y$, for different values of
the channel width. (c) The height direction flow velocity normalized by the
mean flow speed plotted versus channel height (note the rotated orientation of
the graph). For same mean flow speeds, the centerline ($y=w/2$) velocity
varies with the width of the channel. (d) Nondimensional shear rate
($\dot{\gamma}h/\langle u\rangle$) at the wall ($y=0$) plotted versus the
aspect ratio $\mathcal{W}$. For the same mean flow speed $\langle u\rangle$
and height $h$, we can obtain varying shear rates at the wall by changing the
width of the channels.

In our model calculations, $\eta_{c}$ provides a direct measure for the
exclusion zone profile by its definition. From the definition of the
similarity variable, we note that

$\eta_{c}=\frac{Y_{c}}{(3X_{c})^{\frac{1}{3}}}=\frac{y_{c}}{\left(3x_{c}\frac{D_{c}}{\dot{\gamma}}\right)^{\frac{1}{3}}}\leavevmode\nobreak\
,$ (21)

where $(x_{c},y_{c})$ and $(X_{c},Y_{c})$ are, respectively, the dimensional
and nondimensional coordinates of the boundary of the particle excluzion zone.
Therefore, we can write the shape of the boundary in $x,y$ coordinates as

$y=\alpha\eta_{c}(3x(D_{c}/\dot{\gamma}))^{\frac{1}{3}}\leavevmode\nobreak\ ,$
(22)

where $\alpha$ is a prefactor that reflects the details of the three-
dimensional flow.

With the experiment presented in Figure 2, we compare the steady-state
exclusion zone profile with the equation (22) in Figure 5. In section 4, we
used the small $X$ condition
$\displaystyle{\frac{k_{f}}{\dot{\gamma}}X^{\frac{2}{3}}}\ll 1$ to obtain
equation (15). For $k_{f}=0.039$ s-1 and $\dot{\gamma}\approx 2$ s-1 the
condition becomes $x\ll 1.2$ cm, and thus we compare the experimental data
with the equation (22) up to $x=5$ mm. The wall shear rate values were
estimated based on the height-averaged velocity (presented in the next
section) in a channel with the width and height, respectively, $w=1$ mm and
$h=100\leavevmode\nobreak\ \mu$m. In Figure 5, the curves for equation (22)
represent the calculations with wall shear rate values
$\dot{\gamma}=1.74,\leavevmode\nobreak\ 3.47,\leavevmode\nobreak\
8.68,\leavevmode\nobreak\ $ and $17.37$ s-1, corresponding to the mean flow
speeds $\langle u\rangle=50,\leavevmode\nobreak\ 100,\leavevmode\nobreak\
250,\leavevmode\nobreak\ 500\leavevmode\nobreak\ \mu$m/s. The steady-state
exclusion zone match the scaling equation (22) well with the choice of
$\alpha=0.45$. Different channel geometries may require different fitting
parameters. We also obtained the exclusion zone boundary from three-
dimensional (3D) model calculation (details are included in the SI sections 6
and 7); and the scaling analysis, experiments and a 3D model show good
agreement near the inlet of the channel (Figure 5). Note that the 3D
calculation is not the main focus of the study so we only include it in the
SI.

In the experimental situations, the value of $\eta_{c}$ is usually a constant
if the particle zeta potential $\zeta$ and so the mobility $\Gamma_{p}$ are
constant. Therefore, we conclude with the relation that the thickness of
particle exclusion zone $y\propto\dot{\gamma}^{-\frac{1}{3}}$. This means that
the size of the particle exclusion zone is determined by the wall shear rate.
The trend $y\propto\dot{\gamma}^{-\frac{1}{3}}$ was recognized by Lee et al.21
($\delta\propto(U/w)^{-\frac{1}{3}}$; $\delta$ is the exclusion zone thickness
and $U$ is the mean flow speed) in the numerical simulations. In our article,
we further study the rectangular channel flow to calculate the wall shear
rate, and examine the effect of channel geometry and validity of the power law
under multiple conditions.

So far we have confirmed with the experiments and the theoretical model based
on a shear flow approximation that the flow speeds control the exclusion zone
thickness. However, changing flow speed (or the volumetric flow rate) also
changes the number of particles per volume that are excluded from the
boundary. Therefore, in terms of “water cleaning” experiment, changing the
flow rate to test the effect of shear rate on particle exclusion is not a
single-parameter control. Thus in the next section, we study the details of
flow structure in a rectangular microfluidic channel, and show how the wall
shear rate $\dot{\gamma}$ can be varied in such a system without changing the
flow speed. In that way, we can test only the effect of shear rate on the size
of the particle exclusion zone, without perturbing the number of excluded
particles per volume.

##  5 Wall shear rate in a rectangular channel

Fig. 7:  Effect of channel width on the size of the exclusion zone. (a)
Schematic for the width control experiments. Channels with different widths
($w=0.25,\leavevmode\nobreak\ 0.5,\leavevmode\nobreak\
0.75,\leavevmode\nobreak\ 1$ and $2$ mm) were used to demonstrate different
wall shear rates under the same mean flow speed $\langle u\rangle$. (b) Time-
evolution of the exclusion zone measured at $x=8$ mm. (c) Images obtained from
the channels with different widths, at $t=300$ s. For $\langle
u\rangle=100\leavevmode\nobreak\ \mu$m/s and $h=100\leavevmode\nobreak\ \mu$m,
change in $w$ varies the exclusion zone thickness. Horizontal and the vertical
scale bars are, respectively, 1 mm and 500 $\mu$m. (d) Exclusion zone profiles
at $t=300$ s plotted versus $x$. Both (b) time-evolution at a fixed location
and (d) profile at $t=300$ s show larger exclusion zone formation in the
channels with larger width.

In this section, we investigate the flow structure to find important
parameters for controlling the shear rate in a rectangular channel. Consider a
channel with the width, height, and length, respectively $w$, $h$, and $L$. We
define the axes from a corner of the channel (Figure 6(a)), $x$, $y$, and $z$
in the flow, width, and height directions, so that the flow velocity
$\mathbf{u}=u(y,z)\mathbf{e_{x}}$ can be written as

$\displaystyle\begin{split}u(y,z)&=\frac{6\langle
u\rangle}{h^{2}}\left[1-6\left(\frac{h}{w}\right)\sum_{n=0}^{\infty}\lambda_{n}^{-5}\tanh\left(\lambda_{n}\frac{w}{h}\right)\right]^{-1}\\\
&\times\left[(hz-z^{2})-\sum_{n=0}^{\infty}a_{n}\cos\left(\frac{\lambda_{n}\left(z-\frac{h}{2}\right)}{h/2}\right)\cosh\left(\frac{\lambda_{n}\left(y-\frac{w}{2}\right)}{h/2}\right)\right],\end{split}$
(23)

where
$\displaystyle{a_{n}=\frac{h^{2}(-1)^{n}}{\lambda_{n}^{3}\cosh(\lambda_{n}w/h)}}$,
$\displaystyle{\lambda_{n}=\frac{2n+1}{2}\pi}$, and $\langle u\rangle$ is the
mean flow speed. Note that the series solution is first obtained for the
coordinate system defined along the centerline of the channel, then shifted by
$y=w/2$ and $z=h/2$ (see SI for details).

In the experiments, the exclusion zone forms from one side wall, which means
that we are observing the height-averaged view of the system. Therefore, we
consider the height-averaged velocity profile to obtain the wall shear rate at
$y=0$ and simplify the analysis to a two-dimensional flow. The height-averaged
velocity $\bar{u}(y)$ is (see SI for details)

$\bar{u}(y)=\frac{1}{h}\int_{0}^{h}u\leavevmode\nobreak\
dz\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ $ $\displaystyle\begin{split}=\langle
u\rangle&\left[1-\left(\frac{6}{\mathcal{W}}\right)\sum_{n=0}^{\infty}\lambda_{n}^{-5}\tanh\left(\lambda_{n}\mathcal{W}\right)\right]^{-1}\left[1-\sum_{n=0}^{\infty}\frac{6\cosh\left(\frac{\lambda_{n}\left(y-\frac{w}{2}\right)}{h/2}\right)}{\lambda_{n}^{4}\cosh\left(\lambda_{n}\mathcal{W}\right)}\right]\leavevmode\nobreak\
,\end{split}$ (24)

where $\mathcal{W}=w/h$. The equation (24) is normalized by $\langle u\rangle$
and calculated for different values of channel widths ($w=0.25,0.5,0.75$ and
$2$ mm) and $h=100\leavevmode\nobreak\ \mu$m, and plotted versus $y$ in Figure
6(b). Also, the velocity profile viewed from the side wall, along the
centerline, ($u(y=\frac{w}{2},z)/\langle u\rangle$) is calculated and plotted
for channels with different widths (Figure 6(c)). For $w\gg h$, the flow
velocity is almost uniform across the width of the channel, except for the
small region near the walls.

Next we calculate the shear rate $\dot{\gamma}$ at the wall ($y=0$ or $y=w$)
using the height-averaged velocity $\bar{u}$:

$\displaystyle\dot{\gamma}$
$\displaystyle=\frac{d\bar{u}(y)}{dy}\bigg{|}_{y=0,\leavevmode\nobreak\ y=w}$
(25) $\displaystyle=\pm\frac{12\langle
u\rangle}{h}\left[1-\left(\frac{6}{\mathcal{W}}\right)\sum_{n=0}^{\infty}\lambda_{n}^{-5}\tanh\left(\lambda_{n}\mathcal{W}\right)\right]^{-1}\left[\sum_{n=0}^{\infty}\lambda_{n}^{-3}\tanh{(\lambda_{n}\mathcal{W})}\right]\leavevmode\nobreak\
.$

Note that if we maintain the mean flow speed and the channel height constant,
the aspect ratio of the channel $\mathcal{W}$ is the main parameter that can
vary the shear rate at the wall. The normalized shear rate
$\dot{\gamma}h/\langle u\rangle$ is plotted versus $\mathcal{W}$ in Figure
6(d). This calculation shows that for the same mean flow velocity and the
channel height, we can tune the wall shear rate $\dot{\gamma}$ simply by
changing the width of the channel. For smaller $\mathcal{W}$, the shear rate
is higher. We also note that (SI) the shear rates at $z=0$ and $z=h$ are
smaller for larger $\mathcal{W}$, which means that for the same mean flow
velocity and $h$, we obtain lower shear rates on all four walls in the
channels with larger widths.

##  6 Exclusion zone formation in the channels with various geometries

Fig. 8:  (a) Schematic for the height control experiments. Channels with the
width $w=500\leavevmode\nobreak\ \mu$m and various heights were used.
Exclusion zone profiles measured at $t=300$ s in the channels with different
heights, (b) for the same shear rate $\dot{\gamma}$ and (c) for the same mean
flow speed $\langle u\rangle=100\leavevmode\nobreak\ \mu$m/s. (d) Shear rate
$d\bar{u}/dy$ plotted versus $y$ for the channels with different heights.

From the previous sections, we obtained the relation for the thickness of
exclusion zone ($y\sim\dot{\gamma}^{-\frac{1}{3}}$) and calculated the wall
shear rate $\dot{\gamma}$ as a function of channel dimensions. Motivated by
the shear rate calculations, we performed control experiments with the
rectangular channels with the height $h=100\leavevmode\nobreak\ \mu$m and
different widths $w=0.25,\leavevmode\nobreak\ 0.5,\leavevmode\nobreak\
0.75,\leavevmode\nobreak\ 1$ and $2$ mm (Figure 7(a)). We observe the particle
exclusion zone near the CO2-side boundary at a fixed mean flow speed $\langle
u\rangle=100\leavevmode\nobreak\ \mu$m/s. CO2 and N2 pressures were maintained
at, respectively, 10 psi and 1.5 psi. 300 seconds after the CO2 flow is on, we
obtain different exclusion zones in the channels of different widths (Figure
7(b-d)).

As predicted from the calculations, the shear rate is smaller for larger
$\mathcal{W}$ channels with fixed values of $\langle u\rangle$ and $h$, and
thus we observe larger exclusion zones in the channels with larger width. We
measured the time evolution of particle exclusion zone at $x=8$ mm and plotted
versus time in Figure 7(b). The experiments show that the early shape of the
particle exclusion zone also follows the trend predicted by the shear rate
calculations, where we achieve larger exclusion zones in the channels with
larger width. Next we plotted versus $x$ the shape of particle exclusion zone
at $t=300$ s for different channels (distance along the channel) in Figure
7(d). Even considering some overlapping error bars, we observe a clear trend
that the particle exclusion zone is larger in channels of larger width.

We do not directly compare this measurement with equation (22) since the
channel geometries are different and the estimation of $\alpha$ – the
prefactor of equation (22) – based on the experiments done in different
systems may introduce additional errors. Also, the flow structure is
complicated in the wall region, and this complication is quantitatively
different for channels with different dimensions. Therefore, instead, we
further examine the effect of shear rate experimentally by using the channels
with different heights. From the equation (25) we note that the two-
dimensional wall shear rate is a function of $\langle u\rangle$, $h$ and the
aspect ratio $\mathcal{W}$.

In the next set of experiments, we fix the width of the channel $w$ and vary
the height $h$ and the flow speed $\langle u\rangle$. First, experiments were
performed in channels of different heights and with the adjusted flow speeds
to achieve fixed value of $\dot{\gamma}$. For a wide range of the aspect ratio
($1<\mathcal{W}<33.3$), we obtained similar exclusion zones for a single value
of $\dot{\gamma}$. Second, experiments were performed fixing the mean flow
speed $\langle u\rangle$ while varying $h$: we obtain the largest exclusion
zone in the channel with $h=500\leavevmode\nobreak\ \mu$m, and the smallest
for $h=15\leavevmode\nobreak\ \mu$m, which agrees with the trend predicted by
the corresponding shear rates. In Lee et al.21, the excluxion zone was
estimated with the scaling relation $\delta\propto(\langle
u\rangle/w)^{-1/3}$. We show in Figure 8(c) six different systems with the
same values of $\langle u\rangle/w$, which demonstrate that simple estimation
of the shear rate by $\langle u\rangle/w$ may not provide sufficient
information about the two-dimensional exclusion zone forming in channels with
different dimensions. Flow structure and the channel geometry need to be fully
considered.

Fig. 9:  (a,b) Nondimensional, steady-state concentration profiles of CO2 and
ions at $\bar{z}=z/h=0.5$ plotted versus $\bar{x}=x/h$ and $\bar{y}=y/h$.
(c,d) Concentration of ions at $\bar{y}=1$ plotted versus $\bar{z}$ for
different Peclet numbers (note that $\mathcal{W}=10$). When diffusion is
dominant, the ion concentration is more uniform along $\bar{z}$. (e,f)
$\bar{z}$-component of the particle diffusiophoretic velocity plotted versus
$\bar{x}$ and $\bar{y}$. Near the upstream entrance of the CO2 source, highest
$\bar{z}$-direction velocity is created. (g-i) Examining the effect of wall
diffusioosmosis under no channel flow condition. (g) Schematic of the cross-
section of the channel. We solve for the case where $u=0$, and the nonzero
flow velocities $u_{y}$ and $u_{z}$ are induced by wall diffusioosmosis. (h)
Nondimensional, steady-state particle concentration profile plotted versus
$\bar{y}$ and $\bar{z}$. Due to the diffusioosmosis, particle distribution is
not uniform along the $\bar{z}$-direction. (i) Height-averaged $\bar{n}$
plotted versus dimensional $y$, for different channel height values. The
profiles are obtained at $t=300$ s.

We would like to mention more about shear flow approximation in microfluidic
experiments. Figure 6(d) shows $d\bar{u}/dy$ plotted versus $y$. When the
aspect ratio of the channel is large, the change in the height-averaged flow
speed is mostly localized near the wall region, and for small $w/h$, the
change in the shear rate is less dramatic. Therefore, a typical experimental
setup for testing the shear flow approximation is in the orientation where the
observation is made through the small-width window ($\mathcal{W}<1$) by making
tall channels. In that way the experiments can assume a constant shear rate in
the wall region. Knowing this, one can ask about the effect of the orientation
of the channel (whether $\mathcal{W}\geq 1$ or $\mathcal{W}\leq 1$)) on the
water cleaning efficiency, under a fixed volumetric flow rate.

Similar to the height-averaged velocity
$\bar{u}(y)=\frac{1}{h}\int_{0}^{h}u\leavevmode\nobreak\ dz$, we calculated
the width-averaged velocity
$\tilde{u}(z)=\frac{1}{w}\int_{0}^{w}u\leavevmode\nobreak\ dy$ to obtain the
equivalent wall shear rate in the other orientation in the same channel
$\dot{\gamma}_{z}=\frac{d\tilde{u}}{dz}\big{|}_{z=0,\leavevmode\nobreak\ z=h}$
(see SI section 5 for details). When $\dot{\gamma}$ and $\dot{\gamma}_{z}$ are
compared (Figure S5), we obtain that the two-dimensional shear rate is always
larger when the system is viewed from the narrow-side. For example, in a
channel with $w=500$ $\mu$m and $h=100$ $\mu$m, for $\langle
u\rangle=100\leavevmode\nobreak\ \mu$m/s, we obtain $\dot{\gamma}=3.73$ and
$\dot{\gamma}_{z}=6.12$. When the ion concentration gradient is created over
$w$, the measured exclusion zone is expected to be
$(\frac{3.73}{6.12})^{-\frac{1}{3}}\approx 1.2$ times larger than the other
orientation. However, the height is 5 times larger if the ion concentration
gradient is formed across $h$, and thus the volume of particle-free water is
$\approx 4.2$ times larger. An optimization is required for specific design
and fabrication of an actual water cleaning device. Nevertheless, such shear
rate criterion can be a major reference for designing the diffusiophoresis
unit.

So far we have used a shear flow approximation, and, with similarity
transform, we obtained a scaling estimate for the exclusion zone thickness
along the flow. The relation $y\propto\dot{\gamma}^{-\frac{1}{3}}$, along with
the shear rate calculation $\dot{\gamma}=f(\langle u\rangle,h,\mathcal{W})$,
provides useful information for predicting the trend of particle exclusion in
two-dimensional configurations. In the next section, with numerical
calculations we would like to investigate possible three-dimensional effects
that can influence the exclusion zone measurements.

##  7 Three-dimensional effects and wall diffusioosmosis

Consider a three-dimensional rectangular flow (Figure 6(a)) with the flow
velocity defined as equation (23). Then we can solve reaction-diffusion
equations to obtain the CO2 and ion concentrations $c_{c}(x,y,z,t)$ and
$c_{i}(x,y,z,t)$ (details in SI sections 6 and 7).

As a model case, we consider a channel flow with $w=1$ mm,
$h=100\leavevmode\nobreak\ \mu$m, and $\langle
u\rangle=100\leavevmode\nobreak\ \mu$m/s, and plotted the steady-state
concentration profiles of CO2 and ions in Figure 9(a,b) versus $\bar{x}=x/h$
and $\bar{y}=y/x$ (at $\bar{z}=z/h=0.5$). In this configuration, it is often
assumed that the concentration profiles are independent of ${z}$. However, we
note that due to the finite diffusivities of ions, there is a non-uniform
concentration distribution of ions in the ${z}$-direction (Figure 9(c)), near
the upstream of the CO2 source. When the system is diffusion-dominated (Figure
9(d)), the ${z}$-dependence of ion concentration is small. This result
suggests that even though the main concentration gradient is established in
${y}$-direction, there may be a nonzero diffusiophoretic velocity of particles
generated in the height-direction, affecting the particle distribution.

The $\bar{z}$-component of nondimensional diffusiophoretic velocity
$\bar{u}_{pz}=\bar{\Gamma}_{p}\frac{\partial\ln\bar{c}_{i}}{\partial\bar{z}}$
is calculated numerically near the bottom wall ($\bar{z}=0.1$) and plotted
versus $\bar{x}$ and $\bar{y}$ in Figure 9(e,f). Near the inlet of the flow
there are nonzero particle velocities forming in $z$-direction to migrate
particles away from the wall ($\bar{z}=0$). Then for $\bar{x}\gg 1$ this
particle velocity decreases dramatically. A similar trend is observed also in
the width direction; the $\bar{z}$-component of the diffusiophoretic velocity
is highest near the wall where CO2 diffuses in. The nondimensional particle
velocity is a measure for the particle Peclet number by definition
$Pe_{p}=u_{p}h/D_{c}$ and the calculated values suggest that the
$\bar{z}$-direction particle motion is much smaller compared to the main flow
speed ($Pe=5$). Since the ${z}$-direction gradient of ion concentration is
generated due to the channel flow and the finite diffusivity of ions, larger
$\bar{u}_{pz}$ will be obtained for faster flow, and for slow flow
$\bar{u}_{pz}$ will be negligible. This means that the particle exclusion from
the top and bottom walls is expected to be smaller compared to the main
exclusion zone formation in the ${y}$-direction in the presence of the channel
flow. A study of the detailed three-dimensional particle distribution is not
the main focus of the current article and so we include the contour plot of
particle exclusion in the SI. The feature discussed in this section may be of
consideration for maximizing the collection of clean liquid in the in-flow
water cleaning systems.

Fig. 10:  Two-stage water cleaning system. (a,b) Schematics of the channel
geometry and the experiment. (c) Clean water collection in the first stage.
(d) Clean water and particle suspension collection at the end of the second
stage. Scale bars are 1 mm.

One more feature to study in terms of a three-dimensional effect is the
diffusioosmosis on the (top and bottom) walls. When there is a pressure-driven
flow with a moderate speed, analyses that consider particle diffusiophoresis
and the flow velocity may provide sufficient information about the exclusion
zone formation, even though there exists diffusioosmosis forming a slip
boundary condition. We consider one extreme case where the flow velocity is
very small (or zero) and the exclusion zone formation is achieved by mostly
(or only) the diffusion of ions. Small flow velocity decreases the overall
efficiency of the water cleaning system, but if the flow is deliberately
stopped periodically to achieve maximum exclusion zone thickness before each
sample collection step (pulsating flow), one may ask about the influence of
the wall diffusioosmosis.

We consider a model case where there is no flow along the channel and the
migration of particles in the channel is driven by diffusion of ions. We can
solve for the cross-sectional diffusion of ions and diffusiophoresis of
particles and examine the effect of diffusioosmosis (Figure 9(g)). First, we
account for diffusioosmosis along the top and bottom walls, to obtain the flow
velocity ($\mathbf{u}_{f}=u_{y}\mathbf{e}_{y}+u_{z}\mathbf{e}_{z}$) as a
function of $c_{i}$ (more precisely, we use the local chemical equilibrium
approximation 26, 24, and use $c_{c}$ to describe the flow velocity induced by
wall diffusioosmosis; see SI section 8 for details).

In Figure 9(h), the steady-state particle distribution
$\bar{n}(\bar{y},\bar{z})$ is plotted versus $\bar{y}=y/h$ and $\bar{z}=z/h$.
Due to the flow velocity induced by the wall diffusioosmosis, particle
distribution is not uniform along $\bar{z}$ under the concentration gradient
set in the $\bar{y}$-direction. We solved for the particle distribution for
the channel cross-sections with different aspect ratios and plotted the
$\bar{z}$-averaged profiles versus dimensional $y$ (Figure 9(i)). For the
channel with $w=1$ mm, the $\bar{z}$-averaged $\bar{n}$ at $t=300$ s is
plotted for different channel heights. The plots show that at smaller heights,
the effect of diffusioosmosis is larger and exclude particles more. Such
difference in particle exclusion may be significant in microscopic analyses of
diffusiophoretic particle exclusion. However, the slight difference among
different conditions would not affect practical applications, and the
comparison with the no-diffusioosmosis scenario also suggests that the effect
of wall diffusioosmosis can be neglected.

Throughout the paper, we showed by a theoretical model, experiments and
numerical calculations that the thickness of particle exclusion zone near
entrance follows the trend $y\propto\dot{\gamma}^{-\frac{1}{3}}$, where the
shear rate varies with the channel dimensions. By increasing the width and
height of the rectangular channel, we were able to achieve formation of larger
exclusion zone. With these understandings, we would like to propose a possible
two-step design for the in-flow water cleaning system (Figure 10). Again, the
design suggestion is not an optimal geometry or an ultimate solution. We chose
conventional microfluidic fabrication to visualize the concept in a small
length scale. The ideas should be further developed and customized to a real
device application.

##  8 Toward multi-stage water cleaning

In order to visualize the idea of two-stage water cleaning, we fabricated a
PDMS channel ($w_{1}=1$ mm and $h=100\leavevmode\nobreak\ \mu$m) that has the
first bifurcation $\approx 1.5$ cm downstream from the inlet, and the second
bifurcation $\approx 5.5$ cm downstream from the inlet. After the first
bifurcation, the width of the remaining liquid channel is
$w_{2}=750\leavevmode\nobreak\ \mu$m. In order to effectively create the CO2
(or ion) concentration gradients in the second stage, we put a “break” stage
where the liquid channel is in contact with the N2 channel and a channel
filled with the cured optical adhesive (Norland Products; NOA 81 – Figure
10(a,b)). In the break stage, CO2 that is in the liquid phase is expected to
diffuse out on the side facing the N2 channel, as the UV epoxy is impermeable
to gas. While maintaining the direction of the ion concentration gradient, the
absolute concentration of ions decrease at this stage25 (SI). In this way, we
do not perturb the particle distribution by introducing backward ion
concentration gradients. Then in the second stage, an exclusion zone is formed
on the same side as the first stage, and the clean water can be collected at
the final bifurcation outlet. In order to match the size of the exclusion
zones with the model channel design, we adjusted the flow speed to a lower
value ($\langle u\rangle=10\leavevmode\nobreak\ \mu$m/s) than our main
experiments (Figure 10(c,d)).

##  9 Conclusion

Our study on the diffusiophoretically-generated exclusion zone is motivated by
several previous studies that demonstrate diffusiophoresis in microfluidic
systems. By increasing the width of the channels to $O(1)$ mm, we observed
that the exclusion zone forms in a small region near the wall facing the CO2
source. Then, with the shear flow approximation and similarity transform, we
obtained a scaling law for the exclusion zone thickness
($y\propto\dot{\gamma}^{-\frac{1}{3}}$), which led to a detailed calculation
for the shear rate. The wall shear rate in a rectangular channel flow is a
function of $\langle u\rangle$, $\mathcal{W}$ and height $h$, so the channel
geometry and the flow condition can be chosen to obtain a desired wall shear
rate. With control experiments we show that the exclusion zone prediction
using the shear rate is applicable to channels with various widths and
heights. We also confirmed experimentally that a change in the applied CO2
pressure does not significantly vary the exclusion zone thickness in our
system. Finally, we propose a design to demonstrate a continuous, multi-stage
water cleaning. As the most water purification process is never a single step,
by an appropriate choice of the geometrical parameters, such in-flow water
cleaning system can be further scaled up to collect larger amount of cean
water.

## Contributions

S.S. and H.A.S. conceived the project. S.S. designed all experiments with
assist of M.B., and S.S. M.B., and E.H.T. performed the experiments. S.S. and
H.A.S. set up the theoretical model, and S.S. conducted numerical
calculations. M.B. and E.H.T. contributed to the project through the
Independent Work and Princeton Environmental Institute Summer Internship
Program during their undergraduate study at Princeton University.

## Conflicts of interest

S. Shim declares no conflict of interest. H. A. Stone co-founded Phoresis,
Inc., which seeks to utilize CO2-driven diffusiophoresis in the water-cleaning
space.

## Acknowledgements

We acknowledge NSF for support via CBET-1702693. We thank Jesse T. Ault and
Estella Yu for valuable discussions. M.B. and E.H.T. thank the Princeton
Environmental Institute at Princeton University for supporting summer
internships.

## Materials and Methods

The channels are made with PDMS using standard soft lithography. All four
walls are PDMS, and the bottom layer (1-mm PDMS slab) is further attached to
the glass slide for microscopy (Leica DMI4000B). For all experiments using
polystyrene (PS) particles (Thermo Fisher Scientific, diameter = 1 $\mu$m), a
dilute suspension of 0.03 vol$\%$ (in DI water; Milli-Q, EMD Millipore) is
used.

For the experiment with amine-modified polystyrene (a-PS) particles (Fig.
1(d); Sigma Aldrich, diameter 1 $\mu$m), a dilute suspension is prepared at 5
vol$\%$ in DI water. Prior to the experiment, a syringe pump was used to flow
a 1 vol$\%$ aqueous solution of 3-Aminopropyltriethoxysilane (APTES) through
the liquid channel for 15 minutes to prevent the adhesion of amine-modified
polystyrene particles to channel walls 33, 34, followed by flushing with
deionized water for 15 minutes. A N2 gas stream at 1.5 psi is connected for
one minute to dry the central channel.

When running the experiments, the liquid channel was first filled with the
particle suspension, then N2 gas was flowed at 1.5 psi through the top
channel. After one minute of observing steady-state particle flow in the
liquid channel, a CO2 gas stream fixed at 10 psi is connected to the CO2
channel.

The fluorescent images are obtained every 1 second using an inverted
microscope (Leica DMI4000B), with 1.25x magnification. The exclusion zone
boundary was obtained by first dividing the images into 200 separate sectionst
(in $x$-direction), then detecting the mean value of the boundary position in
each section. Each image piece had the width of 6 pixels, and thus the
representative data point is an average of six values. The obtained values
($y$-position of the exclusion zone boundary) are further smoothed by applying
a moving average of period 20. For the time-evolution plots, a similar method
is applied, where the mean value of the boundary position is obtained every
six frames. No smoothing is applied for the time evolution graphs. The smooth
curves in Fig. 1(e) and Fig. S1(b) are obtained using a smoothing spline fit
in Matlab. The image processing and analysis are done with ImageJ and Matlab.

All equations presented in the main text are solved using a combination of
Mathematica and Matlab, and the three-dimensional model included in the SI is
solved with Matlab (see SI for details).

## Notes and references

  * 1 B. V. Derjaguin, G. Sidorenkov, E. Zubashchenko, and E. Kiseleva, Kinetic phenomena in the boundary layers of liquids 1. Capillary osmosis, Kolloidn. Zh., 1947, 9, 335-347.
  * 2 B. V. Derjaguin, S. S. Dukhin, and A. A. Korotkova, Diffusiophoresis in electrolyte solutions and its role in the mechanism of film formation from rubber latexes by the method of ionic deposition, Kolloidn. Zh., 1961, 23, 53-58.
  * 3 S. S. Dukhin, Z. R. Ul’berg, G. L. Dvornichenko, and B. V. Derjaguin, Diffusiophoresis in electrolyte solutions and its application to the formation of surface coatings, Russ. Chem. Bull., 1982, 31, 1535-1544.
  * 4 J. L. Anderson, M. E. Lowell, and D. C. Prieve, Motion of a particle generated by chemical gradients Part 1. Non-electrolytes, J. Fluid. Mech., 1982, 117, 107-121.
  * 5 D. C. Prieve, J. L. Anderson, J. P. Ebel, and M. E. Lowell, Motion of a particle generated by chemical gradients. Part 2. Electrolytes, J. Fluid. Mech., 1984, 148, 247-269.
  * 6 D. C. Prieve, Migration of a colloidal particle in a gradient of electrolyte concentration, Adv. Colloid Interface Sci., 1982, 16, 321-335.
  * 7 D. C. Prieve, and R. Roman, Diffusiophoresis of a rigid sphere through a viscous electrolyte solution, J. Chem. Soc. Faraday Trans. 2, 1987, 83(3), 1287-1306.
  * 8 B. Abécassis, C. Cottin-Bizonne, C. Ybert, A. Ajdari, and L. Bocquet, Nat. Mater., 2008, 7, 785-789.
  * 9 J. J. McDermott, A. Kar, M. Daher, S. Klara, G. Wang, A. Sen, and D. Velegol, Langmuir, 2012, 28, 15491-15497.
  * 10 D. Florea, S. Musa, J. M. R. Huyghe, and H. M. Wyss, Long-range repulsion of colloids driven by ion exchange and diffusiophoresis, Proc. Natl. Acad. Sci. U.S.A., 2014, 111(18), 6554-6559.
  * 11 J. Palacci, C. Cottin-Bizonne, C. Ybert, and L. Bocquet, Osmotic traps for colloids and macromolecules based on logarithmic sensing in salt taxis, Soft Matter, 2012, 8, 980-994.
  * 12 A. Kar, T.-Y. Chiang, I. O. Rivera, A. Sen, and D. Velegol, Enhanced transport into and out of dead-end pores, ACS Nano, 2015, 9, 746-753.
  * 13 S. Shin, E. Um, B. Sabass, J. T. Ault, M. Rahimi, P. B. Warren, and H. A. Stone, Size-dependent control of colloid transport via solute gradients in dead-end channels, Proc. Natl. Acad. Sci., 2016, 113, 257-261.
  * 14 S. Battat, J. T. Ault, S. Shin, S. Khodaparast, and H. A. Stone, Particle entrainment in dead-end pores by diffusiophoresis, Soft Matter, 2019, 15, 3879-3885.
  * 15 N. Shi, R. Nery-Azevedo, A. I. Abdel-Fattah, and T. M. Squires, Diffusiophoretic focusing of suspended colloids, Phys. Rev. Lett., 2016, 117, 258001\.
  * 16 J. T. Ault, P. B. Warren, S. Shin, and H. A. Stone, Diffusiophoresis in one-dimensional solute gradients, Soft Matter, 2017, 13, 9015-9023.
  * 17 A. Gupta, B. Rallabandi, and H. A. Stone, Diffusiophoretic and diffusioosmotic velocities for mixtures of valence-asymmetric electrolytes, Phys. Rev. Fluids, 2019, 4, 043702\.
  * 18 J. L. Wilson, S. Shim, Y. E. Yu, A. Gupta, and H. A. Stone, Diffusiophoresis in multivalent electrolytes, Langmuir, 2020, 36, 7014-7020.
  * 19 A. Gupta, S. Shim, and H. A. Stone, Diffusiophoresis: from dilute to concentrated electrolytes, Soft Matter, 2020, 16, 6975-6984.
  * 20 J. S. Paustian, C. D. Angulo, R. Nery-Azevedo, N. Shi, A. I. Abdel-Fattah, and T. M. Squires, Direct measurements of colloidal solvophoresis under imposed solvent and solute gradients, Langmuir, 2015, 31, 4402-4410.
  * 21 H. Lee, J. Kim, J. Yang, S. W. Seo, and S. J. Kim, Diffusiophoretic exclusion of colloidal particles for continuous water purification. Lab Chip, 2018, 18, 1713\.
  * 22 S. Shin, O. Shardt, P. B. Warren, and H. A. Stone, Membraneless water filtration using CO2, Nature Comm., 2017, 8, 15181\.
  * 23 S. Shim, Interfacial Flows with Heat and Mass Transfer. PhD Thesis, 2017, Princeton University.
  * 24 O. Shardt, B. Rallabandi, S. Shim, S. Shin, P. B. Warren, and H. A. Stone, Driving diffusiophoresis by the dissolution of gas, Under review.
  * 25 S. Shim and H. A. Stone, CO2-leakage-driven diffusiophoresis causes spontaneous accumulation of charged materials in channel flow, Proc. Natl. Acad. Sci. U.S.A., 2020, 117(42), 25985-25990.
  * 26 S. Shim, S. Khodaparast, C.-Y. Lai, J. Yan, J. T. Ault, B. Rallabandi, O. Shardt, and H. A. Stone, CO2-driven diffusiophoresis for maintaining a bacteria-free surface. Soft Matter, 2021, 17, 2568-2576.
  * 27 T. J. Shimokusu, V. G. Maybruck, J. T. Ault, and S. Shin, Colloid separation by CO2-induced diffusiophoresis. Langmuir, 2020, 36, 7032-7038.
  * 28 M. Seo, S. Park, D. Lee, H. Lee, and S. J. Kim, Continuous and spontaneous nanoparticle separation by diffusiophoresis, Lab Chip, 2020, 20, 4118-4127.
  * 29 C. L. Yaws, Handbook of Properties for Environmental and Green Engineering, Gulf Publishing Company, 2008, 1st edn.
  * 30 W. M. Haynes, CRC Handbood of Chemistry and Physics CRC Press, 2015, 96th edn.
  * 31 W. L. Jolly, Modern Inorganic Chemistry McGraw-Hill, New york, 1991, 2nd edn.
  * 32 W. D. Ristenpart and H. A. Stone, Michaelis-Menten kinetics in shear flow: Similarity solutions for multi-step reactions, Biomicrofluidics, 2012, 6, 014108\.
  * 33 Lee, N.Y. and Chung, B.H. Novel Poly(dimethylsiloxane) bonding strategy via room temperature “chemical gluing." Langmuir, 25, 3861–3866 (2009)
  * 34 Zhou, J., Ellis, A.V. and Voelcker, N.H. Recent developments in PDMS surface modification for microfluidic devices. Electrophoresis 31, 2–16 (2010)