Numerical Investigation of a Designed-Inlet Optofluidic Beam Splitter for Split-Angle and Transmission Improvement

The beam splitter is one of the important elements in optical waveguide circuits. To improve the performance of an optofluidic beam splitter, a microchannel including a two-stage main channel with divergent side walls and two pairs of inlet channels is proposed. Besides, the height of the inlets injected with cladding fluid is set to be less than the height of other parts of the microchannel. When we inject calcium chloride solution (cladding fluid) and deionized water (core fluid) into the inlet channels, the gradient refractive index (GRIN) developed in fluids flowing through the microchannel splits the incident light beam into two beams with a larger split angle. Moreover, the designed inlets yield a GRIN distribution which increases the light collected around the middle horizontal line on the objective plane, and so enhances the transmission efficiency of the device. To demonstrate the performance of the proposed beam splitter, we use polydimethylsiloxane to fabricate the microchannel. The results obtained by simulation and experiment are compared to show the effectiveness of the device and the validity of numerical simulation. The influence of the microchannel geometry and the flow rate ratio on the performance of the proposed beam splitter is investigated.


Introduction
Optofluidics aims at the fusion of integrated optics with microfluidics [1]. A wide variety of applications inspired by optofluidics have been developed for solar energy collection [2], water treatment [3], chemical and biological detection and analysis [4,5], sensors [6] and tunable optical devices [7,8]. Optofluidic devices have several advantages over traditional optical devices, including compactness, tunability, reconfigurability and integration. Moreover, the advancements of microfabrication have enabled the realization of some optofluidic devices to be integrated as part of lab-on-a-chip (LOC). Therefore, optofluidics and its applications have earned much attention over the past two decades.
A large number of optofluidic devices have been proposed and investigated, such as microlenses, prisms, waveguides and beam splitters. Seow et al. [9] reported a tunable liquid microlens by controlling the flow rates of the three injected flows with different refractive indices into an expansion chamber. Shi et al. [10] proposed an optofluidic microlens using a static polydimethylsiloxane (PDMS) lens and an air-water interface that can be reshaped to tune the focal point of incident light. Mao et al. [11] reported the design and characterization of a tunable microlens generated by the interface between two miscible fluids of different refractive indices flowing through a 90-degree curve in a microchannel. Mao et al. [12] proposed an optofluidic microlens with liquid gradient refractive index (L-GRIN) to focus light. Mishra et al. [13] presented a liquid-liquid interface microlens formed by two mutually immiscible fluids. By adjusting the hydrostatic pressures and electric field to control the curvature of the interface, the focal length and spherical aberration can be tuned. Le et al. [14] designed an in-plane L-GRIN microlens that realized two-dimensional light beam focusing dynamically. Miccio et al. [15] have shown that a red blood cell can serve as an adaptive optofluidic microlens, which can be used in the application of blood diagnosis. Hamilton et al. [16] developed a three-dimensional focusing device, which was integrated with optical waveguides and used to analyze fluorescent signals from beads in fluid flow. A prism can also be formed by injecting different fluids with different flow rates to a chamber to manipulate light [17,18]. Wolfe et al. proposed an L-GRIN optical waveguide [19] and an optical splitter based on liquid-liquid waveguide [20]. Moreover, with liquids containing different dyes, the device [20] can also be used as a filter. Schmidt and Hawkins presented literature reviews on liquid-core optical waveguides with emphasis on suitability for creating fully planar optofluidic LOC [21] and their fabrication methods and common structures [22]. Yang et al. [23] presented an optofluidic waveguide for light bending and manipulation by pumping two fluids through a straight microchannel, where a bidirectional GRIN-like profile with a higher refractive index in the core region can be achieved at a relatively low flow rate. Yang et al. [24] reported a waveguide with the designed GRIN achieved by the convection-diffusion process of liquid flowing streams with a lower refractive index in the core stream to split a light beam into two parts with a large split angle. It has also been demonstrated the asymmetrical split angle can be achieved by varying the concentration of the cladding flow. Li et al. [25] demonstrated realization of a tunable optofluidic Y-branch waveguides in a planar microchannel by introducing Dean flow.
Waveguides are critical components in optical systems and the optical Y-branch splitter is one of the most important elements in optical waveguide circuits. Yang et al. [24] have shown that a large tunable split angle up to 30 degrees can be achieved by tuning the flow rates of the streams flowing through the inlets of a microchannel based on the designed GRIN achieved by the convection-diffusion process of liquid flowing streams with a lower refractive index in the core stream. In light of Yang's design [24], we try to explore a little further into microchannel structure which yields the GRIN for the optical splitter with a larger split angle and better transmission efficiency. A designed microchannel, including a divergent main channel, two pairs of inlet channels injected with core and cladding fluids and a pair of outlets, is proposed to enlarge the split angle. The L-GRIN beam splitter may enhance the transmission efficiency of the device by setting the height of the inlets injected with cladding fluid to be less than the height of other parts of the microchannel. To demonstrate the performance of the proposed beam splitter, we use PDMS to fabricate the beam splitter. The results obtained by numerical simulation and experiment are compared to show the effectiveness of the proposed device and the validity of numerical simulation. The influence of the geometry parameters of the microchannel and the flow rate ratio of the two injected fluids on the performance of the L-GRIN beam splitter is investigated.

GRIN Beam Splitter Design
To improve the performance of an L-GRIN beam splitter, we adopt a two-stage main channel with symmetric inlet and outlet microchannels and set the sidewalls of each stage of main microchannel to be divergent, as shown in Figure 1. The main channel can be viewed as two parts by taking the centerline of the outlets as the baseline. The primary flow in the first part of the two-stage microchannel is in the direction to the positive y-axis, while that in the second one is in the direction to the negative y-axis. There are two pairs of inlets, A I , B I , A II and B II , and the outlets of the two parts are merged together, as shown in Figure 1a. Here, we use A and B to represent the inlets injected with calcium chloride (CaCl 2 ) solution (cladding fluid) and the inlets injected with deionized water (core fluid), respectively, and use E to represent the outlets. The inlets A I and B I are merged together before entering the main channel; we use C I to represent the merged inlet channel. Similarly, we use C II to represent the merged inlet channel of A II and B II . Thus, the widths of the inlets are denoted by W A I , W B I , W C I , W A II , W B II and W C II , and the width of the outlets is denoted by W E . The subscripts I and II in the symbols used denote the first and the second part of the microchannel, respectively, for convenience. The symbol W I denotes the width of main channel of the first part at the location L I /2 away from the outlet microchannel and we set W I = 2W C I . Similarly, the width of the second part of main channel at the location L II /2 away from the outlet microchannel is denoted by W II and we set W II = 2W C II . To enhance the transmission efficiency of the device, we set the height of the first part of inlets A I and A II , denoted by h, to be less than the height of other parts of the microchannel, denoted by H, and the height of inlets A I and A II is enlarged to H at the locations L S away from the main channel. In Figure 1a, the height of the channels in pink is h and the height of the channels in light blue is H. Cross sections of an inlet branch at the locations with distance from the main channel greater and less than L S are shown in the left and the right of Figure 1b, respectively. In contrast to the bidirectional GRIN achieved by the convection-diffusion process in the beam splitter with uniform channel height proposed by Yang et al. [24], the GRIN distribution achieved by the convection-diffusion process of liquid flowing streams with a lower-refractive-index core stream in the present device with deliberate variation of inlet channel height is three-dimensional (3D). The 3D GRIN, except splits an incident light coupled into the device into two parts, makes more light arrive a zone around the middle horizontal line of an objective plane, where the output optical fibers are coupled to the device. Thus, the 3D GRIN may enhance the transmission efficiency of the present device. The effect of the 3D GRIN generated in the present device on beam splitting will be demonstrated and discussed later.
When we remove the second part of the microchannel, adopt the main microchannel with the side walls parallel to the y-axis, and set the height of inlets (h) to be equal to the height of the main microchannel (H), the present design reduces to a one-stage beam splitter similar to the L-GRIN beam splitter proposed by Yang et al. [24].
There are a very large number of geometrical parameters influencing the performance of the proposed beam splitter. The present work focuses on the effects of geometric modification of the inlets injected with cladding fluid and the two-stage main channel with divergent side walls on the performance of the proposed beam splitter. The former is represented by the ratio of the height of inlets A I and A II to the main-channel height, h/H. The latter is described by the lengths of the two parts of the main channel (L I and L II ) and the divergent angles of the side wall of the two parts of the main channel (β I and β II ).

Numerical Simulations and Geometrical Parameter Analysis
The parameters influencing the performance of the beam splitter considered include h/H, L I , L II , β I , β II and the flow rate ratio of the injected cladding fluid and the injected core fluid. To select the values of those parameters enhancing the performance of beam splitter, we first use numerical simulation to investigate the influence of every single parameter. Then, the efficiency and effectiveness of beam splitters with different combinations of the above parameters are examined. The GRIN distribution formed by the diffusion between the two selected working fluids, calcium chloride solution and deionized water, determines the paths of light in the L-GRIN device. Thus, the simulation includes the convectiondiffusion in the flowing streams and the ray tracing of light beam in the device. Here, we use Computational Fluid Dynamics (CFD) software for the former and self-developed C++ codes for the latter, respectively.
The convection-diffusion process of flowing streams with calcium chloride solution in the cladding stream is described by the continuity equation, the momentum equation and the species convection-diffusion equation for steady flows expressed as where → v is the fluid velocity, in m/s, ρ the fluid density, in kg/m 3 , p the pressure, in N/m 2 , µ the fluid dynamic viscosity, in kg/m · s, D AB the diffusion coefficient of fluid A in fluid B, in m 2 /s, and c A the mass fraction of calcium chloride solution.
In this work, we use the curve-fitting method to construct the relationship between the normalized mass fraction of calcium chloride solution and the fluid density [26]. The result can be expressed as where the normalized mass fraction of calcium chloride solution is defined as with c A,0 denoting the inlet mass fraction of calcium chloride solution. Similarly, the relation between the normalized mass fraction of calcium chloride solution and the fluid viscosity [26] can be expressed as The diffusion coefficient of calcium chloride solution with a mass fraction (c A ) from 0 to 0.3 in deionized water is approximately from 1.153 × 10 −9 to 1.31 × 10 −9 m 2 /s [27]. Here, we take the average value of 1.229 × 10 −9 m 2 /s in the simulation. Calcium chloride solution (fluid A) with c A = c A,0 = 0.3, ρ = 1281.6 kg/m 3 and µ = 0.003467 kg/m · s at 20 • C is injected at inlets A I and A II , and deionized water (fluid B) with c A = 0, ρ = 998.2 kg/m 3 and µ = 0.001003 kg/m · s at 20 • C is injected at inlets B I and B II . The no-slip condition is imposed on all solid walls. The lengths of inlets are set to be ten times of hydraulic diameters based on their cross section and the outlet pressure is set to be atmospheric pressure. The solid walls are set to be impermeable.
In this work, we use the CFD software, ANSYS Fluent 15.0 (ANSYS, Inc., Canonsburg, PA, USA), to simulate the flow field and the concentration distribution by solving the aforementioned governing equations. The governing equations are discretized by using finite volume method. The coupled problems of pressure and velocity in the momentum equation are solved by using Semi-Implicit Method for Pressure Linked Equations (SIM-PLE). The convective terms are handled by using second-order upwind method. During the calculation, the residual at each cell is calculated and normalized by dividing the largest absolute value of the residual in the first five iterations. When the normalized residual is less than 10 −4 , the solution is said to be converged.
After the calculation, we use the CFD Post of ANSYS Fluent to visualize the flow field and export the mass fraction data at every single node which can be converted into the refractive index data. Again, we use curve-fitting method to obtain the relationship between the normalized mass fraction of calcium chloride solution and the refractive index based on the datasheet [26]. The refractive index can be expressed in terms of Φ A (x, y, z) as Thus, we can obtain the refractive index at every grid point, where the value of Φ A is available. In this work, the selected working fluid is the mixture of deionized water (n = 1.33) and calcium chloride solution with a mass fraction of 30% (n = 1.41) and the device is composed by two PDMS layers (n = 1.412). The PDMS has good transparency at visible wavelengths with very low absorption loss. Different L-GRIN distributions can be achieved by adjusting the flow rates of both fluids injected into the microchannel of the optofluidic beam splitter.
When an incident beam from an optical fiber enters the proposed device, it will be refracted continuously by the GRIN medium. To examine the refracted light propagation, we develop ray-tracing codes for the set-up schematically shown in Figure 2. The light source, a green laser with a wavelength of λ= 532 nm, is coupled with an optical fiber which is then set into PDMS layers. The distance between the output end of the optical fiber and the microchannel is L f , and the central lines of the optical fiber and the microchannel are collinear. Moreover, an observation chamber with a vertical wall located at a distance of 60 µm from the side wall of the inlet B II is also fabricated in the set-up, as shown in Figure 2. Here, we use multimode optical fiber as the coupled fiber. The wavelength range of the optical fiber adopted is from 250 to 1100 nm with a refractive index distribution of step index, and the inner and outer diameters of the optical fiber are 100 and 140 µm, respectively. Numerical aperture (N A ) of the optical fiber can be expressed as [28] where n 1 is the refractive index of the core, n 2 is the refractive index of the clad, n is the refractive index of surroundings where the optical fiber is set, and θ A is the half acceptance angle. Since the optical fiber is set in PDMS layers, its numerical aperture is estimated to be 0.16, which leads to a maximum divergent angle of 9 degrees [11]. The greater L f is, the greater the radius of the beam is at the location, where divergent light enters the microchannel. In this work, L f is set to be 50 µm and the radius of the beam at the plane where the rays enter the microchannel is denoted by r 0 . The point of incidence where a ray enters the microchannel on the plane (y = 0) can be written as (x i , 0, z i ) with x i and z i denoting the x and z coordinates of the point of incidence, respectively. To simulate the propagation of the light from an optical fiber, we treat the light energy leaving the outlet of the optical fiber as a bundle belonging to a set of N b bundles of equal energy; the path of the bundle is the same as a ray leaving the optical fiber. On the plane y = 0, the initial position of a ray, (x i , z i ), can be expressed as where R x and R z are random numbers uniformly distributed between 0 and 1, and On account of the divergent angle of the incident beam, we assume that the relationship between incident angle and initial position can be expressed as where θ max is the maximum divergent angle determined by numerical aperture. While assuming geometrical optics holds, the path of rays in the GRIN medium can be obtained by solving the ray equation of geometric optics [29]: where → r = xî + yĵ + zk is the position vector defining a point on the ray, n → r the refractive index at a position → r on the trajectory and s a path length along the trajectory. Nonetheless, the above equation is not in an appropriate form for numerical integration. Employing the variable t defined as dt = ds/n [30], the ray equation can be expressed as Equations (14) and (15) are a system of ordinary differential equations and can be solved with the initial condition that an optical ray vector → T is known and s = 0 at → r = → r i to compute the ray trajectory. Here, we use Runge-Kutta Dormand-Prince (RKDP) method [31] to solve Equations (14) and (15) with a given refractive index distribution. We refer to the series of points obtained in the tracing procedure as the Runge-Kutta points.
While tracing the ray trajectories by the RKDP scheme, refractive-index values and gradients should be given at the Runge-Kutta points. Thus, during the tracing procedure, the refractive-index value and gradient of each Runge-Kutta point have to be retreived from the refractive indices at grid points obtained by the CFD software and Equation (7). Here, we use Moving Least Squares (MLS) method [32,33] to reconstruct the refractive index at the Runge-Kutta points.
When the ray intersects with a surface, the position solved by the RKDP method is usually not exactly at the surface. We use the method based on cubic spline interpolation [34] to determine the ray-surface intersection point. Once we have the coordinates of the ray-surface intersection point, we can calculate the directions of a ray after either reflection or refraction by using the law of reflection or the Snell's law [29].
The description of light propagation shall include the ray trajectories as well as the energy carried by the bundles. When a ray intersects the interface between the working fluid and the PDMS, the reflectivity (ρ) is obtained by Fresnel's relation [29] from the directions of ray and the refractive indices of the media on both sides of the interface. The transmissivity τ at the intersection point equals 1 − ρ. Besides, the energy absorbed by the mixture and the PDMS is assumed to be negligible, and so the energy carried by a bundle only varies when a ray intersects the interface.
To evaluate the performance of the proposed beam splitter, we first calculate the amount of energy carried by bundles arriving the exit plane and the objective planes shown in Figure 2, and then calculate the beam split angle. Here, the amount of energy calculated is normalized with respect to the total energy of bundles leaving the optical fiber. The normalized energy arriving the region located between z = 50 µm and −50 µm on the objective plane is denoted by E ±50µm . Besides, we set a circle with a radius of 50 µm on the objective plane at y = y 1 and calculate the energy gathered inside the circle with its center located at varying x on the positive x-axis. The energy gathered inside the circle is equivalent to the energy arriving an optical fiber coupled with the L-GRIN beam splitter and the radius of the circle is equal to the radius of the optical fiber. Similarly, we normalize the energy gathered inside the circle and denote the normalized energy as E ob (x). When the maximum value of E ob is found in a circle with its center at a location (x 1 , y 1 , 0), the circle is said to be the objective circle and the center of the circle is treated as the center of one of the split beams on the objective plane. Similarly, we can carry out the calculation of the E ob for the circles on the exit plane at y = y 2 to find the center of the split beam, (x 2 , y 2 , 0), on the exit plane. Once we determine (x 1 , y 1 , 0) and (x 2 , y 2 , 0), the split angle can be calculated as

Fabrication and Experiment Apparatus
The microchannel of the proposed beam splitter is symmetric to the middle horizontal plane (z = 0). Thus, the device is constructed by bonding two identical haves of the microchannel. Each half of the microchannel is equipped with two dimensions in height, h/2 and H/2. A master mold of SU-8 (50) film (MicroChem, Round Rock, TX, USA) for a half of the microchannel depicted in Figure 1 is fabricated first by using a two-step lithography process. Next, a PDMS solution (Sylgard 184 elastomer base) and its curing agent are mixed in the weight ratio of 10:1, and then the mixture is poured into the master mold. The workpiece is pumped under vacuum to degas the PDMS and then heated at 100 • C for 60 min in a vacuum drying oven to solidify the PDMS. Subsequently, the PDMS layer is detached from the master mold; each detached PDMS layer is a half of the microchannel. The lower part of the microchannel is stuck on a microscope slide for alignment convenience. Then, the lower part on the microscope slide and the upper part are treated with oxygen plasma. Finally, they are aligned and bonded to finish the fabrication of the whole microchannel. Moreover, a guide channel for an optical fiber delivering light to the microchannel and an observation chamber are also fabricated on the PDMS layers forming the microchannel, as shown in Figure 2.
After the microchannel is fabricated, an F-AM-FC bare fiber adaptor (Newport) is used to connect the light source and an optical fiber. Then, the optical fiber is inserted to the guide channel inside the PDMS layers.
To observe the beam splitting from the top of the microchannel by a microscope connected with a camera, we prepare dyed liquids to be injected into the microchannel by mixing a small amount of Rhodamine B with either calcium chloride solution or deionized water. The fluorescent agent in the dyed liquids serves to visualize the light paths in the microchannel and the observation chamber.
With the silicone tubes attached to the microchannel, we can connect the system to the syringe pumps. To manipulate the L-GRIN in the beam splitter, four streams with different flow rates are required. Thus, we use four syringe pumps to inject calcium chloride solution and deionized water. After the flow achieves steady, we observe the behavior of the flow and the light paths in the L-GRIN beam splitter.

Test of Parameters of Numerical Simulation
To examine the mesh-size independence, we consider a one-stage beam splitter with L I = 240 µm, β I = 0 • and h/H = 0.5. In this work, we are mainly interested in a mixing flow of calcium chloride solution and deionized water yielding higher variation of refractive index. It is found that the diffusion process dominates in the mixing flow with a small flow rate, such as . Q A,I + . Q B,I = 0.067 µL/s. Besides, the diffusion-dominated mixing of working fluids in an optofluidic device leads to an inhomogeneous optical medium with higher variation of refractive index [24]. Thus, we set the total flow rate at 0.067 µL/s and consider the flow rate ratios of the two selected working fluids,  Before applying the ray tracing codes, we select the number of sampled datum points (N data ) and the basis vectors employed by the MLS method for the reconstruction of the refractive indices at the Runge-Kutta points and the number of the incident bundles (N b ) used in the ray tracing procedure for the simulation of light transport. The testing basis vectors include a quadratic, trivariate polynomial expressed as → b 10 → r = 1, x, y, z, xy, yz, zx, x 2 , y 2 , z 2 (18) and a cubic, trivariate polynomial expressed as → b 20 → r = 1, x, y, z, xy, yz, zx, x 2 , y 2 , z 2 , xyz, x 2 y, x 2 z, y 2 x, y 2 z, z 2 x, z 2 y, x 3 , y 3 , z 3 . (19) We consider the light transport through a two-stage microchannel. Based on the comparison of normalized energy fluxes along x-axis with an interval of 5 µm and between z = 2.5 µm and z = −2.5 µm on the objective plane obtained by numerical simulation, we select N data = 64 and → b 10 → r in following simulation. The details of comparison are not presented here for the sake of conciseness. Figure 4 shows the distributions of normalized energy flux along x-direction with an interval of 5 µm and between z = 2.5 µm and z = −2.5 µm on the objective plane obtained by using → b 10 → r , N data = 64, N b = 10 5 , 10 6 and 10 7 . The CPU times of the simulation using N b = 10 5 , 10 6 and 10 7 on i7-4770 CPU @3.40 GHz are 565 s, 5675 s and 56,913 s, respectively. From Figure 4, we can see that the distributions of normalized energy fluxes obtained by using N b = 10 6 and 10 7 are almost overlapping, while the distribution of normalized energy flux obtained by using N b = 10 5 is a little different from the others. Thus, we adopt N b = 10 6 in the cases considered in further discussion.

Comparison of Numerical Results and Experiment Results
To show the effectiveness of the proposed beam splitter and to check the validity of numerical simulation, we fabricate and observe the light transport in a two-stage microchannel with h/H = 0.5, L I = 240 µm, L II = 500 µm, β I = 0 • , β II = 12 • , W A I = 80 µm, W B I = 40 µm, W A II = 150 µm and W B II = 100 µm. Figure 5a shows the top view of ray trajectories in the microchannel filled with deionized water. We can see that the light from source are divergent due to numerical aperture and the rays do not change direction in the medium with uniform refractive index. Then, for comparison purpose, calcium chloride solution is injected at inlets A I and A II and deionized water is injected at inlets B I and B II , respectively, with  Figure 6. We can see that the light paths in Figure 5b and those in Figure 6 are in good agreement. Thus, we use the above numerical parameters in following simulation.

Influence of the Geometry of the Beam Splitter
While investigating the influence of the parameters including h/H, L I , L II , β I , β II and the flow rate ratio of the injected cladding fluid and the injected core fluid, we set W I = 240 µm, W II = 500 µm and H = 240 µm. Thus, the value of h/H corresponds to the height of the inlets injected with cladding fluid (h). In the following simulation, other channel widths are found from the relations, W A I h = W B I H, 2(W A I + W B I = W I , W A II h = W B II H, 2(W A II + W B II = W II , W E = 2W B I for a one-stage microchannel and W E = 2W B I + 2W B II for a two-stage microchannel. To examine the influence of h/H on the GRIN distribution and the light arriving the objective plane, we examine the concentration distribution and light transport in the mixing fluids through a two-stage beam splitter with β I = 0 • , L I = 240 µm, Q B,II = 0.067 µL/s. The designed inlets of the proposed device yield a 3D concentration distribution in the mixing flow through the main microchannel; the gradient of concentration in either the vertical direction (that is, the z direction) or the transverse direction (that is, the x direction) can be seen from the cross-section concentration distributions shown in Figure 7a,b. Equation (7) shows that the refractive index of the mixing streams is directly proportional to the concentration. Thus, the GRIN distributions are corresponding to the concentration distributions for the cases considered; consequently, the refractive index of the liquid core flow stream is lower than that of the liquid cladding flow stream and there is more liquid with higher refractive index around the middle horizontal line (that is, z = 0) of an x-z cross section, as shown in Figure 7. Since light prefers to propagate at a lower phase velocity, the transverse distribution of refractive index makes the input light to be separated and transported from the lower refractive index core region (higher phase velocity) to the higher refractive index cladding region (lower phase velocity), as shown in Figures 5b and 6, and the vertical distribution of refractive index of the mixing streams in the microchannel with an appropriate value of h/H makes more light to arrive the higher refractive index cladding region around the middle horizontal line. The effect of the vertical distribution of refractive index can be confirmed by further simulation of light propagation in the mixing streams flowing in the proposed microchannel. Figure 8 shows the distributions of bundles arriving the objective plane for the cases with various values of h/H. It is found that the distributions of bundles arriving the objective plane for the cases with h/H = 0.5 or 0.4 are preferable to splitting the beam coupled into the central flow stream from an optical fiber. Thus, we adopt the height ratio of 0.5 in following discussion.  Q A,I = 0.067 µL/s. Since the incident rays are divergent and bending to both sides of y-axis in the device, the rays tend to intersect the side walls of main channel with β I < 0 • . Thus, the cases with β I < 0 • are not considered in this work.
The ray trajectories of the cases with various combinations of L I and β I projected on the middle horizontal plane are shown in Figure 9 and the E ob , the E ±50µm and the θ s of those cases are listed in Table 1. From Figure 9, it is readily to see that more rays intersect with the sidewalls of a channel with β I = 0 • , because the rays diverge with propagation. Thus, energy arriving the objective plane decreases with the increase of L I for the cases with β I = 0 • , as shown in Table 1. On the other hand, the E ±50µm and the E ob increases with the increase of L I for the cases with β I = 6 • and 9 • , because the divergent sidewalls of the main channel and the GRIN distribution in flowing streams yield less intersection of the sidewalls and the rays and make more energy bundles arrive the zone between z = 50 µm and −50 µm on the objective plane. From Table 1, we find the following trends. (i) The influence of L I and β I on the E ob is stronger than that of L I and β I on the θ s . (ii) The average of E ob over β I increases with the increase of L I . (iii) The influence of β I on the E ob increases with the increase of L I . (iv) For the cases with L I = 300, 360 or 420 µm, the values of E ±50µm and E ob increase with the increase of β I . (v) In general, the E ob of a case with relatively larger L I and β I is larger. Therefore, we select L I = 420 µm and β I = 9 • as the working parameters of the first part of the microchannel for further parameter study, because using a microchannel with L I = 420 µm and β I = 9 • can yield a high E ob and a medium θ s . Next, based on the analysis of the previous paragraph, we investigate the effects of L II and β II on the performance of beam splitting in the two-stage beam splitter with L I = 420 µm and β I = 9 • . Since more rays entering the second part of the two-stage microchannel have larger bending angles, we consider main channels with L II = 250, 375, 500 µm, β II = 6 • , 9 • , 12 • , and set the flow rate ratio (  Figure 10, we can see that the number of the rays intersecting the channel sidewalls increase with the increase of L II for the cases with a smaller β II (6 • or 9 • ). The intersection of the rays and the sidewalls reduces the light energy arriving the objective plane. Thus, when β II = 6 • or 9 • , the E ±50µm and E ob decrease with the increase of L II , as shown in Table 2. Moreover, when β II = 12 • , both E ±50µm and E ob of the case with L II = 375 µm are maximum among the three cases with β II = 12 • . From Table 2, we also see that the θ s of this case is medium. Similar to the trend found from the flow in the one-stage beam splitter, the influence of L II and β II on E ob is stronger than that of L II and β II on θ s , as shown in Table 2. Therefore, we select L II = 375 µm and β II = 12 o as the working parameters of the second part of the two-stage beam splitter. Using the two-stage beam splitter with L II = 375 µm and β II = 12 • enlarges both of E ob and θ s .

Effect of the Ratio of Flow Rates
In addition to geometry parameters, the ratio of flow rates is the other important parameter, which can be tuned to enhance the performance of the proposed beam splitter. Here, we first investigate the effects of the flow rate ratio in the one-stage beam splitter considered in Section 5.3. Numerical simulation is carried out to calculate the E ob and θ s of the cases with  Table 3, we can see that the increase of the flow rate ratio leads to the decrease of θ s and the increase of E ob , until E ob reaches a maximum value as  Table 4. The concentration distributions on the horizontal plane (z = 0) and the vertical cross sections at y = y 1 = W C I + L I + W E , y = y 1 + L II /2 and y = y 1 + L II of the main channel are shown in Figure 11. From Equation (7), we know that the concentration distribution corresponds to the GRIN distribution. From Table 4, we can see that the increase of . Q A,II / . Q B,II leads to the decrease of θ s . This is because the diffusion zones between the core flow stream and the cladding flow stream get close and the gradient of refractive index decreases with the increase of . Q A,II / . Q B,II , as shown in Figure 11. Such a GRIN distribution reduces the bend angle of the light refracted into the two cladding flow streams. Moreover, from Table 4

Conclusions
An L-GRIN beam splitter with a designed microchannel is proposed to enlarge the split angle and to reduce the light loss of the device. The designed microchannel includes a two-stage main channel with divergent side walls and two pairs of inlet channels, and the height of its inlets injected with cladding fluid is set to be less than the main-channel height. The results of simulation and experiment show the effectiveness of the beam splitter. Some concluding remarks can be drawn from the results obtained.