Similar to the case of solar system planets, reflected starlight from exoplanets is expected to be polarized due to atmospheric scattering, and the net disk-integrated polarization should be nonzero owing to the asymmetrical illumination of the planetary disk. The computation of the disk-integrated reflected flux and its state of polarization involves techniques for the calculation of the local reflection matrices, as well as the numerical recipes for integration over the planetary disks. In this paper, we present a novel approach to calculate the azimuth-dependent reflected intensity vectors at each location on the planetary disk divided into grids. We achieve this by solving the vector radiative transfer equations that describe linear polarization. Our calculations incorporate self-consistent atmospheric models of exoplanets over a wide range of equilibrium temperature, surface gravity, atmospheric composition, and cloud structure. A comparison of the flux and the amount of polarization calculated by considering both single and multiple scattering exhibits the effect of depolarization due to multiple scattering of light depending on the scattering albedo of the atmosphere. We have benchmarked our basic calculations against some of the existing models. We have also presented our models for the hot Jupiter HD 189733 b, indicating the level of precision required by future observations to detect the polarization of this planet in the optical and near-infrared wavelength region. The generic nature and the accuracy offered by our models make them an effective tool for modeling future observations of the polarized light reflected from exoplanets.
