This paper deals with the 3D incompressible nematic liquid crystal flows with density-dependent viscosity in bounded domain. The global well-posedness of strong solutions are established in the vacuum cases, provided the assumption that̄+ ||∇d 0 || L 3 is suitably small with large velocity. Furthermore, the exponential decay of the solution is also obtained. KEYWORDS exponential decay, incompressible nematic liquid crystal flow, strong solutions, vacuum MSC CLASSIFICATION 76A15; 35B65; 35Q35Math Meth Appl Sci. 2020;43:5985-6010. wileyonlinelibrary.com/journal/mma © 2020 John Wiley & Sons, Ltd. 5985 LIUfor some positive constant ,̄. Without loss of generality, both and are normalized to 1. The symbol ∇d ⊙ ∇d, which exhibits the property of the anisotropy of the material, denotes the n × n matrix whose (i, )-th entry is given by i d · d, for i, = 1, 2, 3. We consider an initial boundary value problem for (1) with the following initial and boundary conditions:with d 0 ∈ 2 being given with compatibility, div u 0 = 0 in Ω, and d 0 ∈ C 1 (Ω) satisfying ∇d 0 = 0 on the boundary Ω (see e.g., Hu and Wang 1 and Li 2 ). The nematic liquid crystal flows are regarded as slow-moving particles where the fluid velocity and the alignment of the particles affect each other. The hydrodynamics of nematic liquid crystals is developed by Ericksen 3 and Leslie 4 in the 1960s, but it still retains most important mathematical structures as well as most of the essential difficulties of the original Ericksen-Leslie model. Mathematically, System (1) is a strongly coupled system between the nonhomogeneous incompressible Navier-Stokes equations and the transported heat flows of harmonic map, and thus, its mathematical analysis is full of challenges.When d is a constant vector and |d| = 1, System (1) reduces to the well-known nonhomogeneous incompressible Navier-Stokes equations. Lions 5 , Chapter 2 established the global existence of weak solutions to nonhomogeneous Navier-Stokes equations in any space dimensions for the initial density allowing vacuum. Later, Desjardins 6 proved the global weak solution with more regularity for the two-dimensional case provided that ( 0 ) is a small perturbation of a positive constant in L ∞ norm. Cho-Kim 7 established the local strong solution by using the compatibility condition −div (2 ( 0 )∇u 0 ) + ∇p 0 = √ 0 g and div u 0 = 0 in Ω,