\section{Introduction}
When polymer chains are end-tethered to a surface, a polymer brush is formed, which provides a convenient means to alter the properties of the surface.\cite{mansky97} Introducing nanoparticles (NPs) into the brush gives rise to a wide range of new hybrid materials.\cite{ferhan16,nie16} Such brush-NP systems have broad applications in fields including sensors, information storage devices, medical diagnostic tools, and catalysts, depending on the combination of the NPs and the polymer brush.\cite{christau15} The brush-NP hybrids also serve as a model system to study drug delivery.\cite{mullner16} Furthermore, if the polymer brush is responsive to some perturbations or signals, then the hybrid materials can be used as active components to make various ``smart'' systems.\cite{tokarev12}

For a given brush-NP composite, the state of the brush and the distribution of NPs in the brush are critical factors that determine its prospect for particular applications. These considerations have motivated many experimental studies on controlling the NP distribution and organization in a polymer brush so that certain functions can be realized.\cite{christau15} The hybrid structure is determined collectively by the strength of enthalpic NP-polymer interaction,\cite{diamanti09} the brush characteristics including grafting density,\cite{ferhan12} thickness,\cite{liu02} pattern,\cite{onses12,onses13,steinbach13} composition,\cite{oren09} and possible gradients of these properties,\cite{bhat03,bhat04,bhat06} and the solvent being used.\cite{kesal16}

Various theoretical and computational techniques have been employed to study the interactions between NPs and a polymer brush, including the self-consistent field theory (SCFT),\cite{subramanian96,solis96,steels00,kim02,kim06,kim08,chen05,devos09,devos10,halperin11,egorov12,lian15} Monte Carlo simulations,\cite{milchev08,ermilov10,chen11} molecular dynamics (MD) simulations,\cite{yaneva09,merlitz12,yigit17} dissipative particle dynamics (DPD),\cite{chengjianli15} and Brownian dynamics simulations. \cite{zhang13PhysicaA} These and related methods have also been used to examine various strategies of controlling the distribution and organization of NPs in a polymer brush.\cite{cao07,guskova09,milchev10,opferman12,opferman13,curk13,curk14,zhang13SoftMatter, zhang14JPSB,hua16} For example, Cao and Wu explored the idea of using block copolymer brushes to organize NPs into a 3-dimensional superstructure or a 2-dimensional (2D) layer via a density functional theory.\cite{cao07} Guskova \textit{et al.} observed the cluster formation of polymer-insoluble NPs at the surface of a brush with DPD simulations.\cite{guskova09}  Zhang \textit{et al.} showed via MD simulations that a polymer brush can be used to induce phase separation and crystallization of a binary mixture of NPs embedded in the brush by tuning the size ratio of NPs and the NP-polymer interactions.\cite{zhang13SoftMatter, zhang14JPSB} Hua \textit{et al.} later found that for a similar system the phase separation can also be induced by compressing the NP-containing brush.\cite{hua16}

Despite extensive studies on NP-brush hybrid materials, the effects of solvent evaporation in these systems have been barely explored, even though solvent evaporation is frequently employed in experiments and manufacturing to make polymer-NP composites. Many studies have shown that the evaporation rate plays a critical role in the drying of polymer films and paint.\cite{composto90,strawnhecker01,koombhongse01,luo05,erkselius08,zhang12,kooij15} Jouault \textit{et al.} found that using different casting solvents the same NPs will either aggregate or disperse in the same polymer matrix.\cite{jouault14b} Cheng and Grest showed via MD simulations that the dispersion of NPs in a polymer film is strongly influenced by the evaporation process of the solvent.\cite{cheng16} Kumar \textit{et al.} emphasized the role of solvent evaporation as an outstanding theoretical question in polymer-NP hybrids in a recent review.\cite{kumar2017JCP} Solvent evaporation has also been frequently employed to assemble NPs into various patterns and superstructures.\cite{brinker99,hamon14,josten17,ryu17} It is thus interesting to study the role of solvent evaporation on the structure of brush-NP hybrids.

Here we report on MD simulations of how NPs, initially dispersed in an explicit solvent wetting a polymer brush, are distributed in the brush as the solvent is evaporated. We are particularly interested in determining the condition for a monolayer of NPs to form a well ordered structure (i.e., a hexagonal lattice). One example is shown in Fig.~\ref{evolution}. Our results reveal that depending on brush-NP interactions, various strategies can be used to assemble a well-ordered 2D lattice of NPs, either in or on top of the brush, after solvent evaporation. For strong brush-NP interactions, a polymer brush with a thickness suitable to adsorb a single layer of NPs is needed to yield an ordered 2D array of NPs embedded in the brush. For weak brush-NP interactions, NPs form an ordered 2D array on top of the brush as long as the grafting density is high enough (Fig.~\ref{evolution}). We also find an interesting healing effect where a polymer brush with a low grafting density can cover the grafting plane in a more uniform manner in the presence of weakly adsorbed NPs.

\begin{figure}[ht]
\centering
\includegraphics[width=3.25in]{Fig1.eps}
\caption{Snapshots at various stages during solvent evaporation for the system with $\phi=0.08\sigma^{-2}$, $C=0.65$, and $A_{\rm np} = 100\epsilon$ (see Model and Methodology for the definition of $\phi$, $C$, and $A_{\rm np}$). Color codes: solvent (blue), polymer brush (green), and NPs (orange).}
\label{evolution}
\end{figure}

\section{Model and Methodology}
We model the solvent as a Lennard-Jones (LJ) liquid and the polymer as linear chains of 100 beads connected by permanent bonds. All solvent and polymer beads have a mass $m$ and interact through a standard LJ 12-6 pairwise potential, $U_{\rm LJ}(r)=4\epsilon [(\sigma/r)^{12}-(\sigma/r)^6 -(\sigma/r_c)^{12}+(\sigma/r_c)^6 ]$, where $r$ is the distance between the centers of two beads, $\epsilon$ the strength of interaction, and $\sigma$ the size of beads. The cut off distance is $r_c=3.0\sigma$ for all non-bonded pairs and $r_c=2^{1/6}\sigma$ for bonded pairs of neighboring polymer beads on a chain. The latter are connected by a bond given by the finitely extensible nonlinear elastic potential, $U_{B}(r) = -\frac{1}{2} K R_0^2 {\rm ln}[ 1-( r/R_0 )^2 ]$, where $r$ is the separation distance, $R_0=1.5\sigma$, and $K=30\epsilon/\sigma^2$.\cite{kremer90}

\begin{figure}[ht]
\centering
\includegraphics[width=3in]{Fig2.eps}
\caption{Contact angle $\theta_c$ vs. $A_{\rm np}$ for $D=10\sigma$ (circles) and $\infty$ (squares). Lines are guides to the eye.}
\label{angle}
\end{figure}

The NPs are modeled as spheres of a diameter $d=10\sigma$. Each NP is treated as a uniform distribution of LJ point masses $m$ at a density $1.0m/\sigma^3$, resulting in a mass $M=523.6m$ for the NPs. The NP-NP interaction is described by an integrated form of the LJ 12-6 potential, and its strength is characterized by a Hamaker constant $A_{\rm nn} = 39.48\epsilon$.\cite{Everaers03} To avoid flocculation, the interaction between NPs is chosen to be purely repulsive with the cut off set at $10.438\sigma$.\cite{veld09,grest11} The interaction between a NP and a solvent or polymer bead is determined by a similar integrated potential with Hamaker constants $A_{\rm ns}$ and $A_{\rm np}$, respectively; both interactions are truncated at $9\sigma$. We set $A_{\rm ns}=100\epsilon$, which is sufficiently strong to disperse the NPs in the pure solvent without any polymer chains.\cite{cheng12JCP} We compare two systems with different NP-polymer interaction strengths: $A_{\rm np}=100\epsilon$ and $A_{\rm np}=200\epsilon$. The strength determines the contact angle $\theta_c$ of a NP floating at the surface of a polymer melt that consists of the 100-bead chains modeled here. The results for $\theta_c$ vs. $A_{\rm np}$ are included in Fig.~\ref{angle}. For $A_{\rm np}=100\epsilon$, $\theta_c = 79^\circ$, which indicates that the NPs are only partially adsorbed at the surface of the melt; while for $A_{\rm np}=200\epsilon$, $\theta_c = 0$, which indicates that the NPs initially placed at the surface of the melt eventually diffuse into the bulk of the melt. In Fig.~\ref{angle} we also show the contact angle of a drop formed by 2,346 chains of length 100 on a flat surface, with the same integrated potential for the polymer-surface interaction as that for the polymer-NP interaction extrapolated to $D=\infty$, as a function of $A_{\rm np}$. The data clearly show the effect on $\theta_c$ from the curvature of the NP surface; the smaller the NP, the higher the value of $\theta_c$ at a given $A_{\rm np}$ as long as $0<\theta_c<\pi$.

The model adopted here relies on using LJ potentials and their integrated forms to describe interactions between particles including polymer beads and is thus appropriate for systems dominated by van der Waals interactions. In many experimental systems, NPs and polymers may be charged in a polar solvent such as water and electrostatic interactions may play a major role in controlling the NP-polymer interactions.\cite{bhat03,bhat04,bhat06} In such systems, the NP-polymer interactions may be tuned by varying salt concentrations and exploiting the screening effect of mobile ions on electrostatic interactions. Such situations are of great interest for future studies. Although limited, studies using nonaqueous NP solutions to make NP-brush hybrids were also reported.\cite{oren09,snaith05}

The simulation cell is a rectangular box of dimensions $L_x \times L_y \times L_z$ with $L_x=L_y=100\sigma$. The liquid/vapor interface is parallel to the $x$-$y$ plane, in which periodic boundary conditions are imposed. Polymer chains are randomly end-grafted to the $x$-$y$ plane at $z=0.5\sigma$ with the minimum separation between the tether points set as $r_{tp}$, yielding a planar brush with a grafting density $\phi\equiv N_c/(L_x L_y)$ where $N_c$ is the number of chains. We vary $N_c$ from 100 to 2400; $\phi$ thus changes from $0.01\sigma^{-2}$ to $0.24\sigma^{-2}$. If we take $\sigma \simeq 0.5~{\rm nm}$, then the range of $\phi$ is roughly from $0.04~{\rm nm}^{-2}$ to $0.96~{\rm nm}^{-2}$. In all the simulations, the tether points are fixed and we set $r_{tp}=2\sigma$ except for $N_c=2400$ where $r_{tp}=1.6\sigma$. In the $z$ direction, the system is confined between two flat walls at $z=0$ and $z=L_z$, respectively. The top wall is far above the liquid/vapor interface below which the solvent, NPs, and polymer brush are located. The particle-wall interactions are represented with a LJ 9-3 potential $U_{\rm W}(h)=\epsilon_W [(2/15)(D/h)^{9}-(D/h)^3 -(2/15)(D/h_c)^{9}+(D/h_c)^3 ]$, where $\epsilon_W$ is the interaction strength, $D$ the characteristic length, $h$ the separation between the center of a particle and the wall, and $h_c$ the cut off separation. For all particles, $\epsilon_W=2.0\epsilon$. For the solvent and polymer beads, $D=1.0\sigma$ while for the NPs, $D=5.0\sigma$. We make all the walls repulsive for particles by setting $h_c=0.8583D$ except for the solvent beads at the lower wall, for which $h_c=3.0\sigma$ so that in the initial state the solvent does not dewet the lower wall and the brush.

The systems start in a state where $N$ NPs are randomly dispersed in the solvent that wets the polymer brush and is in equilibrium with the vapor phase of the solvent beads. We define the coverage of the NPs in the $x$-$y$ plane as $C\equiv N\sqrt{3}d^2/(2L_x L_y)$, in which case $C=1$ corresponds a 2D hexagonal close packed layer of NPs in contact that fully covers the $x$-$y$ plane. For the systems studied here, a full coverage with a close packed monolayer requires about 115 NPs. Since the NPs disperse in structures that are not closed packed, we focus primarily on coverages with $20 \le N \le 100$, which correspond to $0.17 \le C \le 0.87$ and can yield full monolayers. The number of solvent beads is varied based on $N_c$ and $N$ to ensure that all the NPs are initially dispersed in the brush-solvent mixture with a volume fraction lower than $10\%$. The approximate location along the $z$-direction for the liquid/vapor interface is estimated and the value of $L_z$ is chosen to give a vapor of thickness of at least $150\sigma$. The systems are then equilibrated for at least $2\times 10^4\tau$, with $\tau \equiv\sqrt{m\sigma^2/\epsilon}$. During equilibration, a vapor phase of the solvent beads is established, a liquid/vapor interface is formed, and the densities of solvent beads in the liquid and vapor phase reach their equilibrium values, $0.62 m/\sigma^{3}$ and $0.055 m/\sigma^{3}$, respectively.

We performed all the simulations using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).\cite{plimpton95,lammps} The equations of motion are integrated using a velocity-Verlet algorithm with a time step $\delta t =0.005\tau$. During the equilibration, the temperature $T$ is held at $1.0\epsilon/k_{\rm B}$ by weakly coupling all beads to a Langevin thermostat with a damping constant $0.1\tau^{-1}$. Once the liquid/vapor interface is equilibrated, the Langevin thermostat is removed except for solvent and polymer beads within $5\sigma$ of the lower wall.\cite{cheng11}

When the evaporation process of the solvent is initiated, the upper wall is moved to $L_z + 20\sigma$ with the region $[L_z, L_z + 20 \sigma]$ designated as the deletion zone. Unless noted, every $0.5\tau$ all the solvent beads in the deletion zone are removed, effectively mimicking the situation that the system is in contact with a vacuum. The evaporation rate which is defined as the number of beads removed every $\tau$ is high initially ($\sim$ $200/\tau$) and decreases with time quickly, approaching a plateau value ($\sim$ $20/\tau$) before diminishing to almost 0 in the late stage of evaporation when most of the solvent is evaporated.\cite{cheng11} We also conducted simulations with fixed, lower initial evaporation rates ($20/\tau$ or $10/\tau$) and found that the structure of the final brush-NP hybrid after solvent evaporation is insensitive to the evaporation rate, at least for the range of evaporation rates accessible to our MD simulations.

\section{Results and Discussion}

\begin{figure*}[ht]
\centering
\includegraphics[width=4.5in]{Fig3.eps}
\caption{Organization of NPs in brushes at various grafting densities after almost all the solvent has evaporated. For all the systems shown, $C=0.65$. For the top three rows (A-C), $A_{\rm np} = 100\epsilon$ while for the bottom three rows (D-F), $A_{\rm np} = 200\epsilon$. (A) and (D) side view; (B) and (E) top view; (C) and (F) Voronoi construction with gray scale based on the value of the Nelson-Halperin order parameter [see Eq.~(\ref{NHEq})] for the NP that each cell encloses. The value of $\phi$ is 0.01, 0.04, 0.08, 0.16, and 0.24$\sigma^{-2}$ for columns 1-5, respectively. In (3E)-(5E), to better show the embedded NPs the polymer beads covering them are not included. The color codes for rows (A), (B), (D), and (E) are the same as in Fig.~\ref{evolution}.}
\label{assembly}
\end{figure*}

Figure~\ref{evolution} illustrates an exemplary process of NP ordering at the surface of a polymer brush during solvent evaporation for $A_{\rm np} = 100\epsilon$. Initially, the NPs are randomly distributed in the solvent, which wets the brush. As the solvent evaporates quickly, the NPs start to accumulate at the receding liquid/vapor interface, because the diffusion of the NPs is much slower than the motion of the interface. Depending on the NP concentration, one or more layers of NPs are formed. For the right NP concentration and a NP-brush interaction leading to $\theta_c > 0$ and favoring phase separation such as in Fig.~\ref{evolution}, one NP layer is formed. In the final state, the NPs are adsorbed at the brush surface and form a 2D triangular lattice with defects.

Simulations such as the one in Fig.~\ref{evolution} are carried out for different combinations of the grafting density $\phi$, the NP coverage $C$, and the strength of the NP-polymer interaction set by $A_{\rm np}$. In Fig.~\ref{assembly} we focus on the effect of varying $\phi$ from 0.01$\sigma^{-2}$ to 0.24$\sigma^{-2}$, and contrast the cases with $A_{\rm np} = 100\epsilon$ ($\theta_c = 79^\circ$) to those with $A_{\rm np} = 200\epsilon$ ($\theta_c = 0$) at a fixed NP coverage $C=0.65$. Embedding a NP in the brush forces the polymer chains around the NP to displace and stretch, which leads to an entropy penalty. For $A_{\rm np} = 100\epsilon$, the enthalpic gain from embedding the NP is not sufficient to compensate for the entropy loss. As a result, in the equilibrium state before evaporation all the NPs are dispersed in the solvent, and the NPs are adsorbed at the surface of the polymer brush after the solvent is evaporated. When the grafting density is low ($\phi = 0.01\sigma^{-2}$), the brush forms islands after drying and the NPs straddle the surface of those islands, as shown in Figs.~\ref{assembly}(1A-1C). As the grafting density is increased to 0.04$\sigma^{-2}$, the brush starts to cover more of the substrate though holes still occur on the surface. The NP distribution after solvent evaporation follows the brush coverage, as shown in Figs.~\ref{assembly}(2A-2C). When the grafting density is further increased to 0.08$\sigma^{-2}$, the brush fully covers the substrate after solvent evaporation and forms a layer with a uniform thickness. After drying the NPs remain on top of the brush and form a lattice in which hexagonal order is easy to identify, as shown in Figs.~\ref{assembly}(3A-3C). Defects still occur and are clearly shown in the Voronoi construction in Fig.~\ref{assembly}(3C). When the grafting density is further increased to 0.16$\sigma^{-2}$ and 0.24$\sigma^{-2}$, as in Figs.~\ref{assembly}(4A-4C) and Figs.~\ref{assembly}(5A-5C), respectively, the situation is similar with the NPs forming a single hexagonal layer on top of the brush after the solvent is evaporated. The lattice usually still contains defects which can move around or disappear over time. However, the removal of a defect typically requires long time, indicating that defects can be quite persistent.

Systems with $A_{\rm np} = 200\epsilon$ behave differently before and during solvent evaporation. For $A_{\rm np} = 200\epsilon$ the NP-brush interaction is strong enough to compete with the entropic penalty associated with embedding a NP in the brush. In the equilibrium state before evaporation, the number of NPs engulfed in the polymer brush grows as $\phi$ increases for the range of $\phi$ studied here and is expected to decrease as $\phi$ increases further, reflecting a balance between the enthalpic gain from NP-polymer attractions and the entropic penalty of embedding NPs in the brush. At a low grafting density $\phi = 0.01\sigma^{-2}$, the NPs are adsorbed onto the domains formed by polymer chains and tend to penetrate the regions with a thickness comparable to the NP diameter after solvent evaporation. In the regions where the brush thickness is larger than the NP diameter, some NPs are embedded in the brush and can even form more than one layers, as shown in Fig.~\ref{assembly}(1D). When the grafting density is quadrupled to 0.04$\sigma^{-2}$, the coverage of brush on the substrate increases but holes still form where there are neither polymer chains nor NPs. In the area coated by the brush, most NPs are engulfed by polymer chains and form a layer close to the substrate [Fig.~\ref{assembly}(2D)]. Since not all the NPs can be accommodated by the brush which only partially covers the substrate, some NPs are adsorbed on top of the first layer and usually stay in registry with the first layer (i.e., the projection of those NPs falls between the NPs in the first layer). The Voronoi constructions shown in Figs.~\ref{assembly}(1F) and (2F) are for the first layer of NPs. Although the overall packing quality is low, some local hexagonal order can still be observed.

The situation is quite interesting at a grafting density $\phi = 0.08\sigma^{-2}$, where the brush is dense enough to fully cover the substrate and its thickness is similar to the NP diameter [Fig.~\ref{assembly}(3D)]. In this case the NPs are all fully embedded in the brush after the evaporation of the solvent and form a single layer with almost perfect hexagonal order as shown in Figs.~\ref{assembly}(3E-3F). Matching the brush thickness and the NP coverage is important for obtaining the well ordered monolayer as can be seen in the data for thicker brushes. When the grafting density is further increased, the polymer chains in the brush are more stretched because of the steric repulsion between the chains, which overcomes the entropy penalty of chain stretching. As a result, the brush thickness is increased to more than twice the NP diameter, which means that the brush is able to engulf more than one layers of NPs. Figs.~\ref{assembly}(4D-4F) and Figs.~\ref{assembly}(5D-5F) show results for $\phi = 0.16\sigma^{-2}$ and $0.24\sigma^{-2}$, respectively. At $C=0.65$, some NPs penetrate deep into the brush while more NPs are embedded in a region close to the brush surface. These latter NPs form a layer in which hexagonal packing is frequently seen but defects are abundant, as shown in the Voronoi constructions for the top layers in Figs.~\ref{assembly}(4F) and (5F). The packing order in this layer is much worse than that in the system with $\phi = 0.08\sigma^{-2}$ [Fig.~\ref{assembly}(3F)]. If the NP concentration in the solution before evaporation is increased to $C \gg 1$, then a 3-dimensional lattice is expected for NPs embedded in a brush with a thickness much larger than the NP diameter. Such a scenario is interesting for future studies.

Solvent evaporation is intrinsically an out-of-equilibrium process and as a result the final distribution of NPs in the polymer brush observed in our simulations may not be an equilibrated one. In particular, our data show that in the final stage of drying, NPs diffuse very slowly inside the brush for strong NP-brush interactions and dense brushes. Instead, only local vibrations of NPs around their average locations are observed in simulations. The results in Figs.~\ref{assembly}(4D-4F) and Figs.~\ref{assembly}(5D-5F) may thus actually correspond to kinetically trapped states. However, these states do not seem to evolve much on time scales accessible to MD simulations and even in additional annealing runs.

\begin{figure}[ht]
\centering
\includegraphics[width=3in]{Fig4.eps}
\caption{The fraction of NPs having 6 neighbors, $f_6$ (open symbols), and the Nelson-Halperin order parameter, $m_6$ (solid symbols), (a) as a function of $\phi$ when $C$ is fixed at 0.65, and (b) as a function of $C$ when $\phi$ is fixed at $0.08\sigma^{-2}$. The data are for the systems with $A_{\rm np} = 100\epsilon$ (circles) and $A_{\rm np} = 200\epsilon$ (squares), respectively. Lines are guides to the eye.}
\label{bond_order}
\end{figure}

To quantify the packing order, i.e., how close a lattice is to an ideal hexagonal lattice, we have computed two parameters. One is $f_6$, the fraction of NPs having exactly 6 neighbors. Another is the Nelson-Halperin order parameter for bond orientation defined as\cite{nelson79}
\begin{equation}
m_6 =  \left| \frac{1}{N_l}\sum_{j=1}^{N_l}
\frac{1}{n_j}\sum_{k=1}^{n_j}\exp (i6\theta_{jk})\right| ,
\label{NHEq}
\end{equation}
where the first sum is over all NPs in the lattice, the second sum is over all the $n_j$ nearest neighbors of the $j$-th NP, and $\theta_{jk}$ is the angle formed by the bond between the $j$-th NP and its $k$-th nearest neighbor with respect to a fixed axis. From the definition $f_6=1$ and $m_6=1$ indicate an ideal hexagonal lattice and smaller values correspond to lattices with disorder and defects. Results for $f_6$ and $m_6$ for the systems in Fig.~\ref{assembly} are shown in Fig.~\ref{bond_order}(a). The two parameters follow the same trends as $\phi$ varies. In particular, for $A_{\rm np} = 200\epsilon$ both $f_6$ and $m_6$ peak at $\phi=0.08\sigma^{-2}$, indicating that a polymer brush with an optimized grafting density is required to produce a well-packed 2D hexagonal lattice of NPs [see Fig.~\ref{assembly}(3D)] when the NP-brush interaction is strong. For $A_{\rm np} = 100\epsilon$, however, the value of $f_6$ and $m_6$ are close to 1 as long as $\phi\gtrsim 0.08\sigma^{-2}$, indicating that a dense brush is capable to adsorb a layer of NPs on its surface and the NPs naturally order into a hexagonal lattice with few defects.

In Fig.~\ref{packing} the effect of NP coverage on the resulting lattice of NPs after solvent evaporation is shown for $\phi=0.08\sigma^{-2}$. The top two rows are for $A_{\rm np} = 100\epsilon$. At a low NP coverage ($C=0.43$), the packing is random and no hexagonal order is observed as shown in Fig.~\ref{packing}(1B). When the NP coverage is increased to $C=0.65$, a hexagonal lattice emerges after the evaporation of the solvent. Note that the results in Figs.~\ref{packing}(2A-2B) are obtained in a second, independent run from the one used to generate Figs.~\ref{assembly}(3A-3C). The fact that in all of our simulations the partially adsorbed NPs order into a hexagonal packing on top of the brush after evaporating the solvent shows that the result is quite robust. When the NP coverage is further increased to $C=0.87$, some NPs diffuse into the brush but the rest form a hexagonal lattice with few defects at the surface of the brush, as shown in Figs.~\ref{packing}(3A-3B), confirming that the hexagonal order occurs at the brush surface as long as the NP coverage is high enough and the brush is dense enough to yield a relatively uniform thickness. The results on $f_6$ and $m_6$ in Fig.~\ref{bond_order}(b) show the same trend.

\begin{figure}[ht]
\centering
\includegraphics[width=3.25in]{Fig5.eps}
\caption{The organization of NPs at various coverages on top of a brush with $\phi=0.08\sigma^{-2}$ after solvent evaporation. (A) and (C) top view; (B) and (D) Voronoi construction. (1) $C=0.43$, (2) $C=0.65$, and (3) $C=0.87$. For the three systems in (A) and (B), $A_{\rm np} = 100\epsilon$. For (C) and (D), $A_{\rm np} = 200\epsilon$. Note that in (1C)-(3C) all NPs are embedded in the brush; for better visibility a layer of polymer beads at the top is removed. Note that (2C) and (2D) here are identical to Fig.~\ref{assembly}(3E) and Fig.~\ref{assembly}(3F), respectively. They are duplicated to make comparison more convenient. The color codes are the same as in Fig.~\ref{assembly}.}
\label{packing}
\end{figure}

\begin{figure*}[ht]
\centering
\includegraphics[width=6in]{Fig6.eps}
\caption{The adsorption of NPs at various coverages in a brush with $\phi=0.04\sigma^{-2}$ after solvent evaporation. (A) side view and (B) top view. (1) dry brush, (2) $C=0.17$, (3) $C=0.43$, (4) $C=0.65$, (5) $C=0.87$, and (6) $C=0.43$. The NP-polymer interaction strength is $A_{\rm np} = 100\epsilon$ for (2)-(5) while $A_{\rm np} = 200\epsilon$ for (6). The color codes are the same as in Fig.~\ref{assembly}.}
\label{healing}
\end{figure*}

In Fig.~\ref{bond_order}(b) and the bottom two rows of Fig.~\ref{packing}, the results for the systems with $\phi=0.08\sigma^{-2}$ and $A_{\rm np} = 200\epsilon$ are also included. For $C=0.43$, all the NPs are embedded at similar heights in the brush, forming a single layer. However, the NP coverage is low and only local hexagonal packing occurs within disorder at large scales as shown in Fig.~\ref{packing}(1C). For $C=0.65$, a hexagonal NP lattice with few if any defects is formed. For $C=0.87$, about 80\% of the NPs penetrate into the brush and form a layer close to the substrate. The Voronoi construction for this first layer is shown in Fig.~\ref{packing}(3D). The remaining NPs are randomly distributed above the first layer. Therefore, for $A_{\rm np} = 200\epsilon$ a 2D hexagonal lattice of NPs is only obtained when both $\phi$ and $C$ are fine tuned to yield a system where a single layer of NPs with high enough loading is embedded in the brush. The same trend is also shown in the quantitative results in Fig.~\ref{bond_order} for $f_6$ and $m_6$ as functions of $\phi$ and $C$, respectively. A brush with a lower $\phi$ than the optimal value cannot fully cover the whole substrate [e.g., Figs.~\ref{assembly}(2D-2F)] while at a higher $\phi$ the brush is thick enough to support two layers of NPs [e.g., Figs.~\ref{assembly}(5D-5F)] or more. A system with $C$ lower than optimal has NPs forming a single layer in the brush but the layer is too dilute for hexagonal packing [e.g., Figs.~\ref{packing}(1C-1D)]. For systems with $C$ higher than optimal, NPs form more than one layers with diminished or no hexagonal order in each layer [e.g., Figs.~\ref{packing}(3C-3D)]. It should be pointed out that when $C$ is higher than optimal, the distribution of NPs engulfed in the brush varies with $\phi$. In Figs.~\ref{packing}(3C-3D) for $\phi=0.08\sigma^{-2}$ and $C=0.87$, the majority of the NPs first forms a layer close to the substrate with the rest distributed above this layer. However, when $\phi$ is increased to $0.16\sigma^{-2}$ and $0.24\sigma^{-2}$, most of the NPs engulfed in the brush form a layer close to the brush surface with the rest distributed below this layer, as shown in Figs.~\ref{assembly}(4D-4F) and Figs.~\ref{assembly}(5D-5F) for $C=0.65$. The same phenomenon is also observed for $\phi=0.24\sigma^{-2}$ and $C=0.87$ with a major layer of NPs near the brush surface (results not shown). This change is consistent with the idea that there is a larger entropic penalty of penetrating a denser polymer brush.

Finally we have observed an interesting effect brought about by the adsorption of NPs on a polymer brush with a low grafting density. The results in Fig.~\ref{healing} are for a brush with $\phi = 0.04\sigma^{-2}$ at various values of NP coverage. For the dry brush (i.e., $C=0$), the chains are not sufficient to cover the whole substrate and form domains that percolate the plane of the substrate as shown in Figs.~\ref{healing}(1A-1B); many holes can be seen. At $C=0.17$, the brush covers most of the surface with only one hole. When the NPs are adsorbed on top of the brush for the case of a weak NP-polymer interaction (e.g., $A_{\rm np} = 100\epsilon$), the thickness becomes more uniform. When the NP coverage is increased to $C=0.43$, the hole completely disappears and the brush is able to fully cover the substrate, which we call a ``healing'' effect of the adsorbed NPs. Interestingly, the hole reappears when the NP coverage is increased to larger values [Figs.~\ref{healing}(4A-4B) and (5A-5B)], indicating the more NPs do not strengthen the ``healing'' effect, though the brush is partially ``healed'' if we compare Figs.~\ref{healing}(4B) and (5B) to Fig.~\ref{healing}(1B). One possible issue is that the high evaporation rates resulting from evaporating the solvent into a vacuum may play a role in the ``healing'' capability of partially adsorbed NPs. Within our limited range of evaporation rates, we find the same behavior, but much slower rates would be an interesting direction for future studies.

The ``healing'' effect was never observed for strong NP-polymer interactions. For example, for $A_{\rm np} = 200\epsilon$ ($\theta_c = 0$) the brush is still highly nonuniform for $C=0.43$ [Figs.~\ref{healing}(6A-6B)], which shows a good ``healing'' capability as shown in Figs.~\ref{healing}(3A-3B) for $A_{\rm np} = 100\epsilon$ ($\theta_c = 79^\circ$). These results reveal an interesting strategy to make a low-grafting-density brush more uniform by using adsorbed, but not embedded, NPs at intermediate values of coverage.

\section{Conclusions}
In this paper we have presented results from large-scale molecular dynamics simulations on ordering NPs with a polymer brush via solvent evaporation. It has been found that a dense brush with grafting density $\phi \gtrsim 0.08\sigma^{-2}$ and a weak NP-polymer interaction resulting in NP adsorption on top of the brush can be combined to robustly yield a hexagonally packed lattice of NPs after solvent evaporation as long as the NP coverage $C$ is in the range leading to a single layer that is not too dilute. If the NP-polymer interactions are strong, then a system with an appropriate combination of $\phi$ and $C$ is required to yield a single layer of NPs engulfed in the brush with a thickness similar to the NP diameter. This layer of NPs has a well-ordered hexagonal arrangement. Our results thus point to strategies to generate a 2D hexagonal lattice of NPs using a polymer brush depending on the interactions between the NPs and the brush. Our results also reveal an interesting effect that the adsorption of NPs can ``heal'' a polymer brush with a low grafting density and make the brush to be more uniform and even fully cover the substrate. This discovery may lead to interesting applications of brush-NP hybrid systems as a protective layer of a solid surface.

\section*{Acknowledgments}
Acknowledgment is made to the Donors of the American Chemical Society Petroleum Research Fund (PRF \#56103-DNI6), for support of this research. This research used resources of the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the United States Department of Energy under Contract No. DE-AC02-05CH11231. This work was performed, in part, at the Center for Integrated Nanotechnologies, an Office of Science User Facility operated for the U.S. Department of Energy (DOE) Office of Science. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA-0003525.


\section{Introduction}
The propagation of electrons has many similarities with the propagation of light.\cite{Datta,W.van}
In two-dimensional electron gas (2DEG), where the mean free path is larger
than the size of the system, ballistic electrons propagate following straight-line trajectories
which is analogous to light rays.
When the ballistic electrons transmit across a semiconductor junction, electrons should undergo
refraction in analogy with light refraction at optical boundary with different refractive indices.\cite{U.Sivan,H.van,J.Spector,Lv}
Such phenomena can be understood simply in terms of Snell's law, where the refractive index of photons is replaced by the wave vector
of electrons.
Thus, it is possible to manipulate electrons like photons and electron optics has attracted worldwide attention because of its underlying applications.
Electron focusing, diffraction, and double-slit interference experiments are examples of electron optics which have been clearly observed in 2DEG systems.\cite{add1,add2}

In conventional 2DEG systems, electrostatic lenses have been demonstrated in high mobility GaAs about thirty years ago.\cite{U.Sivan,J.Spector}
Since then, many work have been undertaken to obtain various electron optical devices like mirrors, prisms, lenses and splitters.\cite{A1,A2,A3,A4,A5}
One interesting topic of electron optics is the negative refraction,\cite{B1,B2}
which is challenging to achieve in conventional optical systems. For photons, this behavior can be realized in optical metamaterials.\cite{B3,add3}
In electronic systems, negative refraction can be achieved quite straightforward.\cite{A1,B4}
For example, when electrons transmit across a p-n junction, in order to conserve the transverse component of momentum, the transverse group velocity has to change a sign between the valence bands in p side and the conduction bands in n side, hence leads to the negative refraction.
Graphene has a unique band structure with a linear dispersion relation of low-lying
excitations, which gives rise to many peculiar properties.\cite{c3,ref1} Because of the vanishing band gap and the high intrinsic mobility, graphene has also been considered as an attractive platform for studying the electron optics.\cite{A1,A2,A3,A4,c1,c2,c4,ref2}
Recent experiments have clearly demonstrated the negative refraction in graphene p-n junctions.\cite{B4,c4}
This negative refraction can be used to design a perfect lens, and has many other potential applications.

In optics, a beam of light at the anisotropic crystal interface exhibits the double
refraction effects. Recently, monolayer transition-metal dichalcogenides (TMDs) have been successfully fabricated in experiments.\cite{IsSC1,IsSC2,IsSC3,IsSC4,IsSCth}
The TMDs (e.g. NbSe$_2$, MoS$_2$) exhibit Ising pairing in superconducting phase at sufficiently low temperature, with an in-plane
upper critical field far above the Pauli paramagnetic limit.\cite{IsSC1,IsSC2,IsSC3,IsSC4} This unusual behavior is attributed to the inter-valley pairing protected
by Zeeman-type spin-valley locking against external magnetic fields.\cite{IsSC2} This out-of-plane Zeeman-type spin polarization of the valleys can
be used to achieve double refraction in electron optics. Motivated by this,
in this paper, we propose a model to realize double refraction and double focusing
of electric current by a normal-hexagonal semiconductor junction.
This model is based on a hexagonal lattice system (like graphene or TMDs), but breaks the A-B lattice symmetry.
By introducing a Rashba spin-orbit interaction into the system,
the two valleys in the Brillouin zone are no longer equivalent.
When an electron is incident from the normal conductor side with the square lattice,
double refraction and double focusing occur, and the electron is transmitted to two different modes.
The two modes have different refractive indices due to the two valleys being inequivalent.
We investigate electron optics in both n-n junction and p-n junction below.
In particular, the incident electron undergoes double negative refraction in p-n junction.
We also show that by double negative refraction with elaborately tuning the chemical potential of the system,
the p-n junction can work as a spin splitter (valley splitter), which may be useful in spintronics
(valleytronics).

The rest of this paper is organized as follows.
In Sec. \ref{sec:Hamiltonian}, we present the model Hamiltonian and the corresponding band structures of the normal-hexagonal semiconductor junction, considering both the effect of
A-B lattice symmetry breaking and Rashba spin-orbit interaction.
In Sec. \ref{sec:Double refraction}, we investigate the positive double refraction in n-n junction and the negative double refraction in p-n junction.
In Sec. \ref{sec:Spin splitter}, a spin splitter is designed based on the negative double refraction.
Finally, a brief summary is presented in Sec \ref{sec:conclusions}.

\begin{figure}[!htb]
\includegraphics[width=1.0\columnwidth]{1.eps}
\caption{(Color online) Lattice diagram and band structures.
(\textbf{a}) The schematic diagram of the normal-hexagonal semiconductor junction in the lattice model.
(\textbf{b}) Brillouin zone of the hexagonal lattice.
(\textbf{c}) Energy bands along the dash line in (\textbf{b}) in which $k_y=0$.
Dashed red lines denote the energy bands of spin-up electrons, while the blue lines are the spin-down bands.
(\textbf{d}) Zoom-in figure of the energy bands of the dotted box in (\textbf{c}).
The parameters are $\epsilon_{R}=-6.2$, $t_{R} =-8$, $\beta=4$, and $\beta_{s}=-0.03$.}
\label{fig1}
\end{figure}

\section{\label{sec:Hamiltonian}Model Hamiltonian and band structures}
In the following, we introduce the model Hamiltonian and demonstrate the corresponding band structures.
Fig. \ref{fig1}a is the lattice model of the normal-hexagonal semiconductor junction,
which consists of a square lattice at the left side and a hexagonal lattice at the right side.
The Hamiltonian of the whole system can be written as
\begin{align}\label{eq:a}
H=H_{L}+H_{R}+H_{T},
\end{align}
where $H_{L}$, $H_{R}$, $H_{T}$ are the Hamiltonians of the normal, hexagonal semiconductors,
and the coupling between them, respectively.
For the normal conductor, we consider the square lattice model with the dispersion relation of its carrier being quadratic. In the tight-binding representation, the Hamiltonian $H_{L}$ is of the form\cite{d1,d2}
\begin{align}\label{eq:b}
H_{L}=\sum_{i\sigma}\epsilon_{L} a^{\dagger}_{i\sigma} a_{i\sigma}+ \sum_{\langle{ij}\rangle \sigma} t_{L} a^{\dagger}_{i\sigma}a_{j\sigma},
\end{align}
where $a_{i\sigma}$ and $a^{\dagger}_{i\sigma}$ are the annihilation and creation operators
at the discrete site $i$, and $\epsilon_{L}$ is the on-site energy.
The second term in equation (\ref{eq:b}) is the nearest-neighbor hopping. $t_{L}$ is
the hopping energy, which is positive for valence bands and negative for conduction bands.
Experimentally, it is possible to break the A-B lattice symmetry in a hexagonal lattice like graphene.
For example, isolated graphene/BN bilayers break the chemical equivalence of graphene A and B lattice sites.\cite{E1,E2} Graphene growth on the reconstructed surface of MgO(111) also leads to A-B lattice symmetry breaking.\cite{E3,E4}
In addition, spin-orbit interaction can also play a very important role in some two-dimensional materials.
The emergence of the quantum spin Hall effect and topological insulators can be attributed to
the rise of the spin-orbit interaction.\cite{F1,F2,F3,F4}
Thus, we consider a general Hamiltonian at the hexagonal lattice side as
\begin{align}\label{eq:c}
H_{R}=\sum_{i\sigma}(\epsilon_{R}+\lambda_{i}\beta) b^{\dagger}_{i\sigma} b_{i\sigma}+ \sum_{\langle{ij}\rangle\sigma} t_{R} b^{\dagger}_{i\sigma}b_{j\sigma}\nonumber \\
+\sum_{\langle\langle{ij}\rangle\rangle \sigma \sigma^{'}}i\beta_{s}\nu_{ij}s_{\sigma \sigma^{'}}^{z}b_{i\sigma}^{\dagger}b_{j\sigma^{'}},
\end{align}
where $b_{i\sigma}$ and $b^{\dagger}_{i\sigma}$ are the annihilation and creation operators
at the discrete site $i$ of the right side. $\epsilon_{R}$ is the on-site energy,
and $\beta$ represents the energy difference between A-B sublattice.
Here $\lambda_{i}=\pm1$ for A(B)-sublattice.
The second term is the nearest-neighbor hopping term, and $t_{R}$ is the hopping energy.
The third term is the spin-orbit interaction which connects second nearest neighbor.
The same term also appears in the seminal work of the quantum spin Hall effect in graphene.\cite{F1}
$s_{z}$ is a Pauli matrix representing the electron's spin and $\nu_{ij}=-\nu_{ji}=\pm1$ depending on the orientation of the two site $i$ to $j$.\cite{F1}
The spin-orbit coupling $\beta_{s}$ is usually very small comparing to $t_{R}$,
but it can be in the order of meV by the Bi-cluster deposition.\cite{ref3}
The Hamiltonian $H_{T}$ of the coupling between the left and right lead is
\begin{align}\label{eq:d}
H_{T}=\sum_{ij\sigma}t_{c}a_{i\sigma}^{\dagger}b_{j\sigma}+h.c.
\end{align}
where $t_c$ is the coupling strength. Here we assume that the interface of the normal
and hexagonal semiconductors are perfect (see Fig.1a).
In usual, there are defects in the real devices and the interface is rough as well.
While the incident and refractive electrons are in the bottom of the conduction band
or the top of the valence band, their wavelengths are usually very long.
As long as the wavelength is much longer than the defect size and the interface roughness,
the direction of the refraction can be well maintained and the results are barely affected
by the defect and the imperfect interface.

The Brillouin zone of the hexagonal lattice is demonstrated in Fig.\ref{fig1}b.
We choose each unit cell containing four atoms (see red
box in Fig.\ref{fig1}a) which can simplify the calculations of the transmission coefficients,
hence the Brillouin zone is one half smaller than the usual Brillouin zone of graphene.
In Fig.\ref{fig1}c, we plot the energy bands along the dashed line in Fig.\ref{fig1}b with $k_y=0$.
Because of the Brillouin zone being half smaller, the number of the energy bands is double, and each spin orientation has four energy bands.
Breaking the A-B lattice symmetry induces a big energy gap at the K and K' points.
We mainly consider the electrons transport near the Fermi surface,
and the four lower energy bands can be ignored, as they are much below the Fermi surface $E=0$.
In Fig.\ref{fig1}d, the zoom-in figure of the energy bands near the K and K' points are presented.
The spin-orbit coupling can be viewed as an effective magnetic field $\bf{B}$ that points in the opposite directions at the K and K' points.
At K point, spin up electrons have lower energy while at K' point the other way around.
Thus, for energy bands with a definite spin orientation, the energy dispersion relation between the K and K' points are no longer equivalent any more.
In fact, this energy band is the same as that of TMDs in the normal phase.\cite{IsSCth}
So a junction consisting of a normal conductor coupled with the TMDs can be regarded
as a real example of our model. Double refraction occurs when electrons incident from normal conductor with the square lattice side
transmit to the two valleys, as demonstrated in the next section.

\begin{figure}[!htb]
\includegraphics[width=1.0\columnwidth]{2.eps}
\caption{(Color online) Heterojunction band diagrams and schematic of double refraction.
(\textbf{a}) The band diagrams of the n-n junction.
Incident electron i transmits to the two valleys with the
transmission coefficients $T_{1}$ and $T_{2}$.
(\textbf{b}) The corresponding p-n junction case.
(\textbf{c}) and (\textbf{d}) are the corresponding equienergy lines in the $k_x$-$k_y$ plane for
the n-n junction in (\textbf{a}) and the p-n junction in (\textbf{b}).
Here the arrows indicate the directions of the incident and refraction carriers.
(\textbf{e}) Schematic of the double positive refraction in the n-n junction. Electrons emitted
from the source i are bent different amounts at the n-n junction interface, forming two virtual focuses
at the same side of the electron source. Only one virtual focus is plotted in (\textbf{e}).
(\textbf{f}) Schematic of the double negative refraction in the p-n junction.
Electrons emitted from the source i undergo double negative refraction at the p-n junction interface,
forming two real focus at the opposite side of the electron source.}
\label{fig2}
\end{figure}

\section{\label{sec:Double refraction}Double refraction and double focusing}
In this section, we demonstrate how the normal-hexagonal semiconductor junctions leads
to the double refraction and double focusing.
Due to the Pauli matrix $\sigma_z$ commuting with the Hamiltonian in Eq.(\ref{eq:a}),
the spin in the $z$ direction is conserved in the scattering process.
So in the following analysis, we first consider a spin component, e.g. spin-up one.
Fig.\ref{fig2}a and \ref{fig2}b respectively illustrate the band diagrams of
the n-n junction and p-n junction for a spin component,
and the corresponding equienergy lines in the $k_x$-$k_y$ plane are shown in Fig.\ref{fig2}c and \ref{fig2}d.
For an incident carrier (electron or hole) from the left normal conductor, there are two
beams of outgoing electrons, due to the two non-equivalent valleys.
Notice that the directions of the two beams of outgoing electrons are usually different.
So the double refraction occurs in this normal-hexagonal semiconductor junction.
Because of the transverse translation invariance in the model,
the transverse component of the momentum ($k_y$) is conserved when electrons transmit across the junction.
So the transverse group velocity preserves the sign between the n-n junction (see Fig.\ref{fig2}c),
but changes a sign between the valence bands in the p side and the conduction bands in the n side (see Fig.\ref{fig2}d).
As a consequence, the double refraction is positive in a n-n junction,
while negative in a p-n junction.
Fig.\ref{fig2}e and \ref{fig2}f are the schematic diagrams of the positive and
negative double reflections, respectively.
For the positive double refraction, the outgoing electrons
are bent different amounts as they pass through the junction and form two virtual focuses
at the same side of the electron source i.
However, in Fig.\ref{fig2}f, two real focuses are formed at the opposite side of the electron
source, attributing to the opposite transverse velocity and the negative refractive index.

\begin{figure}[!htb]
\includegraphics[width=1.0\columnwidth]{3.eps}
\caption{(Color online) Transmission coefficients and refractive indices for the n-n junction.
(\textbf{a}) and (\textbf{b}) are the transmission coefficients $T_{1}$ and $T_{2}$
at the two valleys as a function of the incident angle $\theta$ for the different incident energy $E$.
(\textbf{c}) and (\textbf{d}) are the corresponding refractive indices.
The parameters are $\epsilon_{L}=7$, $t_{L} =-8$, others are the same as that in Fig.\ref{fig1} .}
\label{fig3}
\end{figure}

To further investigate the nature of the double refraction, we calculate the transmission coefficients and
the refractive indices for both case in detail.
By applying the non-equilibrium Green's function method to the tight-binding Hamiltonian $H$
in equation (\ref{eq:a}), the transmission coefficients $T_1$ and $T_2$
can be obtained quite straightforward by following the same procedure as in Refs.[\onlinecite{G2,G1}].
Fig.\ref{fig3} displays the case of the double positive refraction in a n-n junction,
where (a) and (b) are the transmission coefficients and (c) and (d) are the refractive indices.
For small incident angle, the transmission coefficient $T_{1}$ is nearly a constant,
and is insensitive to the incident energy $E$.
As the incident angle increases, $T_{1}$ rapidly reduces and $T_{2}$ significantly increases.
For large incident angle, both $T_{1}$ and $T_{2}$ drop to zero, attributing to the absence of the
corresponding transmition modes at the n side. Thus, electrons undergo total reflection at large incident angle.
Note that $T_{1}$ dominates the whole scattering process for a large range of parameters.
The increasing of the incident energy enhence $T_{1}$ at large incident angle,
but meanwhile suppress the transmission coefficient $T_{2}$.
The total transmission probability $T=T_{1}+T_{2}\approx 1$ indicates that the n-n junction is
almost transparent for the electrons.
Meanwhile, the two transmission modes have different refractive property as shown in Fig.\ref{fig3}c and \ref{fig3}d.
Here the refractive index is define as $n_{i}=\sin{\theta}/\sin{\theta_{i}}$, where $\theta$ is the incident angle
and $\theta_{i}$ ($i=1,2$) is the refraction angle of transmission mode $i$.
The refractive index $n_{1}$ is insensitive to the incident energy, and increases as the incident angle grows.
We emphasize that $n<1$ for a large range of parameters,
hence the electronic refraction at this valley serves as optically thinner medium in optics.
For a perfect virtual focus, the refractive index n obeys the law $n^2=c+(1-c)\sin^2\theta$,
where the parameter $0<c<1$ in optically thinner medium.
The refractive index $n_{1}$ approximately follows this law, hence this virtual focus is nearly perfect
(see Fig.\ref{fig2}e).
For the refraction at the other valley, the refractive index $n>1$ for a large range of parameters,
hence electronic refraction at this valley serves as optically denser medium in optics.
Here $c>1$ in the optically denser medium. Although $n_{2}$ is sensitive to incident energy,
it also approximate matches the law $n^2=c+(1-c)\sin^2\theta$.
Therefore, the second virtual focus is also relatively good
(The second virtual focus is not plotted in Fig.\ref{fig2}e).

\begin{figure}[!htb]
\includegraphics[width=1.0\columnwidth]{4.eps}
\caption{ (Color online) Transmission coefficients and refractive indices for the p-n junction.
(\textbf{a}) and (\textbf{b}) are the transmission coefficients $T_{1}$ and $T_{2}$ at the two valleys
versus the incident angle $\theta$ for the different incident energy $E$.
(\textbf{c}) and (\textbf{d}) are the corresponding refractive indices.
The parameters are $\epsilon_{L}=-7$, $t_{L} =8$, others are the same as that in Fig.\ref{fig1} .}
\label{fig4}
\end{figure}

Next, we focus on the p-n junction in which the double negative refraction occurs,
as illustrated in Fig.\ref{fig4}.
The main features of the transmission coefficients are basically the same as the case of positive refraction, except the position of $T_{1}$ and $T_{2}$ are changed.
In this case, $T_{2}$ dominates the scattering process (see Fig.\ref{fig4}a and \ref{fig4}b),
but $T_{1}$ is not very small and it has the same order of $T_2$.
As for the refractive indices shown in Fig.\ref{fig4}c and \ref{fig4}d,
$n_{1}$ and $n_{2}$ become negative and almost have a mirror symmetry with the positive refractive
indices in Fig.\ref{fig3}c and \ref{fig3}d about the axis $n=0$.
Contrary to the previous case, two real focus are formed at the opposite side of the electron source,
and the focusing performance of the incident electrons is pretty good (see Fig. \ref{fig2}f).

The spin-down electrons show similar refractive behaviors at the two valleys. The double refraction is positive
for a n-n junction and negative for a p-n junction, and the absolute value of the refractive index $|n|>1$ at one valley and $|n|<1$ at the other.
Because of the time reversal invariance of our Hamiltonian, the scattering matrix of the spin-up and spin-down electrons
satisfy the relation ${\bf s}_{\uparrow}^{T}={\bf s}_{\downarrow}$, which indicates that $T_{1\uparrow}+T_{2\uparrow}=T_{1\downarrow}+T_{2\downarrow}$.

\begin{figure}[!htb]
\includegraphics[width=1.0\columnwidth]{5.eps}
\caption{(Color online) A spin splitter based on the normal-hexagonal p-n junction.
(\textbf{a}) is the band diagrams of the spin splitter based on the normal-hexagonal p-n junction.
Here the Fermi energy $E$ locates between the energy bands of the up-spin and down-spin species
of the n side.
(\textbf{b}) is the schematic diagram of the electron's trajectory in the spin splitter.
Electrons with different spin are separated in real space and focus in different regions. }
\label{fig5}
\end{figure}

\section{\label{sec:Spin splitter}Spin splitter and detective device}

Spin-orbit coupling modifies the electronic band structures at the two valleys and
the double refraction become spin dependent as mentioned above.
Thus, for spin non-polarized incident electrons from the left side,
the outgoing electrons on the right side will be spin-polarized,
i.e. the normal-hexagonal semiconductor junction can be worked as a spin splitter.
Below, as an example, we consider the p-n junction system, and the energy bands are shown in Fig.\ref{fig5}a.
Here the red lines denote the energy bands of spin-up electrons on the n side while the blue lines denote spin-down ones, and the green lines are the bands on the p side where the spin-up and spin-down electrons are degenerate.
In order to clearly illustrate the physical picture of the spin splitter,
the Fermi energy $E$ is set to locate between the red line and blue line (see Fig.\ref{fig5}a).
In this case, the incident spin-up electrons can only get through the junction via valley K at the state $T_{1\uparrow}$
while the spin-down electrons via valley K' at the state $T_{2\downarrow}$.
Notice the states $T_{1\uparrow}$ and $T_{2\downarrow}$ are not mutually the time-reversal states and
the dispersion relations at the two valleys $\epsilon_{K\uparrow}({\bf k}-{\bf K}) \not= \epsilon_{K'\downarrow}({\bf k}-{\bf K'})$,
although the system has the time-reversal invariance with
the dispersion relations having $\epsilon_{K\uparrow}({\bf k}-{\bf K}) = \epsilon_{K'\downarrow}(-{\bf k}-{\bf K'})$.
So the refractive indices of the two valleys are different, leading that the electrons with
different spin are separated in real space and focus in different regions.
By detailed calculation of the incident and refraction angles,
we plot the schematic diagram of the electron's trajectory in Fig.\ref{fig5}b.
Spin-up electrons focus in a small area while spin-down electrons focus in a relatively large area.
The electron's outgoing area can be divided in three regions.
In region 1, both spin-up and spin-down electrons can be observed,
but spin-down electrons are in the majority.
In region 2, only spin-up outgoing electrons could be observed.
In region 3, we can hardly detect any outgoing electron.
In particular, at two focal points the density of the spin-up and spin-down outgoing electrons are very high,
which can be detected by using ferromagnetic STM.
Note that here the spin-up and spin-down electrons are at the valleys K and K', respectively. So
the spin splitter is also a valley splitter .

\begin{figure}[!htb]
\includegraphics[width=1.0\columnwidth]{6.eps}
\caption{(Color online) Transmission coefficients and refractive indices for the spin splitter.
(\textbf{a}) and (\textbf{b}) are the transmission coefficients $T_{1\uparrow}$ and $T_{2\downarrow}$ for spin-up and spin-down electrons, respectively.
(\textbf{c}) and (\textbf{d}) are the corresponding refractive indices.
The parameters are $\epsilon_{R}=-4.2$, $t_{R} =8$, $\beta=4$, $\beta_{s}=-0.1$, $\epsilon_{L}=-7$ and $t_{L} =8$.}
\label{fig6}
\end{figure}

Let us study the transmission coefficients and refractive indices of the spin splitter.
Fig.\ref{fig6}a and \ref{fig6}b are the transmission coefficients $T_{1\uparrow}$ and $T_{2\downarrow}$ of
the two transport channels.
The two transmission coefficients are in fact identical for various incident energy.
This is becasue for each spin there is only one transmission mode, and the general relation
$T_{1\uparrow}+T_{2\uparrow}=T_{1\downarrow}+T_{2\downarrow}$ reduces to $T_{1\uparrow}=T_{2\downarrow}$
in present case.
$T_{1\uparrow}$ and $T_{2\downarrow}$ are the maximum when the incident angle $\theta=0$.
With the increase of $\theta$, they decrease. And at large $\theta$,
$T_{1\uparrow}$ and $T_{2\downarrow}$ drops rapidly to zero.
However, due to the band structures at the two valleys being not completely the same,
the refractive indices present different behavior as illustrated in Fig.\ref{fig6}c and \ref{fig6}d.
With larger incident angle, the refractive index $n_{1}$ decreases, while $n_{2}$ shows opposite behavior.
The refractive index $n_{1}$ approximately follows the asymptotic behaviour of perfect lenses,
and the spin-up electrons focus in a small area.
However, $n_{2}$ obviously do not follow the law, and the spin-down electrons focus in a relatively large area.
Note that by changing the spin-orbit coupling constant $\beta_{s}$,
different types of focusing pattern can be realised.

\section{\label{sec:conclusions}Summary}
In conclusion, we have proposed a model to realize double refraction
and double focusing of electric current using a normal-hexagonal semiconductor junction.
By breaking the A-B lattice symmetry and introducing the Rashba spin-orbit interaction,
a large energy gap is formed and spin-up and spin-down electrons experience opposite effective magnetic field at the two valleys.
Incident electrons transmit to the two valleys, and result in a double refraction at
the hexagonal lattice side.
The double refraction can be either positive or negative, depending on the junction being n-n type or p-n type.
The two valleys reveal different type of refractive behaviors.
The absolute value of the refractive index $|n|>1$ at one valley and $|n|<1$ at the other.
Additionally, the p-n junction can be used as a spin splitter (valley splitter),
which could be verified by the ferromagnetic STM.
Our results may be useful for the engineering of double
lenses in electronic system and have underlying application in spintronics.

\section*{Acknowledgement}
This work was financially supported by National Key R and D Program of China (2017YFA0303301),
NBRP of China (2015CB921102) and NSF-China under Grants Nos. 11574007 and 11274364.

\section*{References}