Ultrasound elastography is a prominent noninvasive medical imaging technique that estimates tissue elastic properties to detect abnormalities in an organ. A common approximation to tissue elastic modulus is tissue strain induced after mechanical stimulation. To compute tissue strain, ultrasound radio frequency (RF) data can be processed using energy-based algorithms. These algorithms suffer from ill-posedness to tackle. A continuity constraint along with the data amplitude similarity is imposed to obtain a unique solution to the time-delay estimation (TDE) problem. Existing energy-based methods exploit the first-order spatial derivative of the displacement field to construct a regularizer. This first-order regularization scheme alone is not fully consistent with the mechanics of tissue deformation while perturbed with an external force. As a consequence, state-of-the-art techniques suffer from two crucial drawbacks. First, the strain map is not sufficiently smooth in uniform tissue regions. Second, the edges of the hard or soft inclusions are not well-defined in the image. Herein, we address these issues by formulating a novel regularizer taking both first- and second-order derivatives of the displacement field into account. The second-order constraint, which is the principal novelty of this work, contributes both to background continuity and edge sharpness by suppressing spurious noisy edges and enhancing strong boundaries. We name the proposed technique: Second-Order Ultrasound eLastography (SOUL). Comparative assessment of qualitative and quantitative results shows that SOUL substantially outperforms three recently developed TDE algorithms called Hybrid, GLUE, and MPWC-Net++. SOUL yields 27.72%, 62.56%, and 81.37% improvements of the signal-to-noise ratio (SNR) and 72.35%, 54.03%, and 65.17% improvements of the contrast-to-noise ratio (CNR) over GLUE with data pertaining to simulation, phantom, and in vivo tissue, respectively. The SOUL code can be downloaded from code.sonography.ai.