Heat Assisted Magnetic Recording (HAMR) is a developing data-storage technology in which a laser delivery system is integrated to the conventional HDD air bearing slider that carries the read and write transducers. The laser beam heats a small spot of around 20nm size on the storage media up to few hundred degrees Celsius [1]. This heating causes several effects on the lubricant such as temperature gradient, thermocapillary shear stress, viscosity drop, and evaporation followed by its depletion. Conventionally [2–8], the disk lubricant is considered as a Newtonian viscous fluid that can be fully described by a viscosity parameter μ. However, in rapid heating and forcing conditions like HAMR, the time dependent nature of the lubricant becomes very important. Measurements [9] show that under some conditions the lubricant behaves like a Maxwell viscoelastic fluid that can be described by two parameters: viscosity μ and Maxwell relaxation time λ. Itoh et al. [10] show that the viscoelastic behavior becomes even more considerable in the case of sub-10nm material confinement. Both the Maxwell relaxation time and viscosity can be functions of temperature and lubricant thickness in case of ultra-thin film lubrication. Karis [9] measured Maxwell relaxation time as a function of temperature for a variety of lubricants, such as Z-dol and Z-tetraol. Fig. 1 shows the results of these measurements for both Z-dol and Z-tetraol. Maxwell relaxation time plays a vital role in determining the behavior of the material under thermal and mechanical loads. In order to have a proper understanding of the effect of Maxwell relaxation time, we non-dimensionalize this parameter by the timescale of the problem to introduce a non-dimensional Deborah number De = λU/L. So, De is a function of HAMR temperature T, disk speed U, laser spot size L, and lubricant type. For purely-viscous materials both the Maxwell relaxation time and De are zero and for purely-elastic materials, both are infinity. So in the case of viscoelasticity, if De ≪ 1 the viscosity mode is dominant, if De ≫ 1 the elasticity is dominant, and if De ≈ 1 the material behaves viscoelastically. Therefore, De is good measure for the viscoelastic behavior of the material. Some attempts have been made to fit the lubrication theory for viscoelastic materials using perturbation methods. However these methods require that the Deborah Number be small enough [11]. Fig. 2 shows the Deborah Number as a function of laser spot size for different lubricant temperatures. Accordingly, at the target of a HAMR laser spot size of L = 20nm, the Deborah number is very large and therefore, the material behaves less viscous and more elastic. Therefore, the traditional methods of lubrication theory cannot describe the lubricant’s behavior in this limit. Consequently, we developed anew direct Finite Element Method (FEM) approach to simulate the behavior of the linear viscoelastic Maxwell fluid lubricants under HAMR conditions.