We investigate higher-order derivative estimates for the Lamé system with hard inclusions embedded in a bounded domain in \(\mathbb {R}^{d}\) . As the distance \(\varepsilon \) between two closely spaced hard inclusions approaches zero, the stress in the narrow regions between the inclusions increases significantly. This stress is captured by the gradient of the solution. The key contribution of this paper is a detailed characterization of this singularity, achieved by deriving higher derivative estimates for solutions to the Lamé system with partially infinite coefficients. These upper bounds are shown to be sharp in two and three dimensions when the domain exhibits certain symmetries. To the best of our knowledge, this is the first work to precisely quantify the singular behavior of higher derivatives in the Lamé system with hard inclusions. We also consider the system with soft inclusions and obtain upper bounds for higher-order derivatives.
