In this paper, we employ the linear virtual element spaces to discretize the semilinear sine-Gordon equation in two dimensions. The salient features of the virtual element method (VEM) are: (a) it does not require explicit form of the shape functions to construct the nonlinear and the bilinear terms, and (b) relaxes the constraint on the mesh topology by allowing the domain to be discretized with general polygons consisting of both convex and concave elements, and (c) easy mesh refinements (hanging nodes and interfaces are allowed ). The nonlinear source term is discretized by employing the product approximation technique and for temporal discretization, the Crank-Nicolson scheme is used. The resulting nonlinear equations are solved using the Newton's method. We derive a priori error estimations in L 2 and H 1 norms. The convergence properties and the accuracy of the virtual element method for the solution of the sine-Gordon equation are demonstrated with academic numerical experiments. studies focused on undamped sine-Gordon equation. The damped sine-Gordon equation demands more sophisticated numerical techniques. In this direction, the following articles are presented in the literature: Nakajima et al. [6] proposed an efficient numerical technique for the solution of the sine-Gordon equation considering dimensionless loss factors and unit less normalized bias. Recently, a dual reciprocity boundary element method [4] and the radial basis function method [7] have been introduced for the above mentioned model problem. However, most of these techniques require the domain to be discretized with simplex elements such as triangles/quadrilaterals.The introduction of elements with arbitrary number of sides has revolutionized the simulation technique. Elements with arbitrary number of sides offer greater flexibility in meshing and are less sensitive to mesh distortion [8,9]. Among the different approaches based on polygonal elements, such as the scaled boundary finite element method [10], the polygonal FEM [8,11,12], the strain smoothing techniques [13,14,15], the virtual element method (VEM) [16,17,18] has an edge. The VEM is a recently developed framework, which provides FEM-like approximations over arbitrary polygonal meshes. Some of the salient features of the VEM are: (a) it does not require an explicit form of the shape functions to compute the bilinear form; (b) it can be applied to convex, and concave elements, i.e., star convexity is not required; (c) easy mesh refinements (hanging nodes and interfaces are allowed ) and (d) the degrees of freedom associated with the VEM space, provides sufficient information to compute the discrete bilinear form. Due to its versatility, the method has been applied to a wide variety of problems. Some of them include: the Stokes equation [19], the Plate bending equation [20], linear elliptic, parabolic, and hyperbolic equations [16,17,18,21], linear and nonlinear elasticity problems [22,23], the convection diffusion equation with small diffusion [24], eigenvalue problems [25,26], ...