In this work, we develop a two-phase lattice Boltzmann method (LBM) to simulate axisymmetric thermocapillary flows. This method simulates the immiscible axisymmetric two-phase flow by an improved color-gradient model, in which the single-phase collision, perturbation and recoloring operators are all presented with the axisymmetric effect taken into account in a simple and computational consistent manner. An additional lattice Boltzmann equation is introduced to describe the evolution of the axisymmetric temperature field, which is coupled to the hydrodynamic equations through an equation of state. This method is first validated by simulations of Rayleigh-Benard convection in a vertical cylinder and thermocapillary migration of a deformable droplet at various Marangoni numbers. It is then used to simulate the thermocapillary migration of two spherical droplets in a constant applied temperature gradient along their line of centers, and the influence of the Marangoni number (Ca), initial distance between droplets (S-0), and the radius ratio of the leading to trailing droplets (A) on the migration process is systematically studied. As Ma increases, the thermal wake behind the leading droplet strengthens, resulting in the transition of the droplet migration from coalescence to non-coalescence; and also, the final distance between droplets increases with Ma for the non-coalescence cases. The variation of,S-0 does not change the final state of the droplets although it has a direct impact on the migration process. In contrast, A can significantly influence the migration process of both droplets and their final state: at low Ma, decreasing A favors the coalescence of both droplets; at high Ma, the two droplets do not coalesce eventually but migrate with the same velocity for the small values of A, and decreasing A leads to a shorter equilibrium time and a faster migration velocity. (C) 2016 Elsevier Ltd. All rights reserved.