A high-order local transmitting boundary is developed to model the propagation of elastic waves in unbounded domains. This transmitting boundary is applicable to scalar and vector waves, to unbounded domains of arbitrary geometry and to anisotropic materials. The formulation is based on a continuedfraction solution of the dynamic-stiffness matrix of an unbounded domain. The coefficient matrices of the continued fraction are determined recursively from the scaled boundary finite element equation in dynamic stiffness. The solution converges rapidly over the whole frequency range as the order of the continued fraction increases. Using the continued-fraction solution and introducing auxiliary variables, a high-order local transmitting boundary is formulated as an equation of motion with symmetric and frequency-independent coefficient matrices. It can be coupled seamlessly with finite elements. Standard procedures in structural dynamics are directly applicable for evaluating the response in the frequency and time domains. Analytical and numerical examples demonstrate the high rate of convergence and efficiency of this high-order local transmitting boundary.Rigorous methods are, in general, spatially and temporally global. They are computationally expensive for large-scale problems and long-time calculations. One of the most popular rigorous methods is the boundary element method [11,12]. This method satisfies the governing equations in the problem domain and the radiation condition at infinity automatically by using a fundamental solution. Only the boundary needs to be discretized. However, it is very complicated to evaluate the fundamental solution when the material is anisotropic. The thin-layer method (consistent boundary condition) [13] has been developed for horizontally layered media. In this method, the boundary in the vertical direction is discretized and displacement functions in the horizontal direction satisfy the radiation condition at infinity. Exact non-reflecting boundary conditions, such as Dirichlet-to-Neumann (DtN) map [14], can be constructed from analytical solutions for unbounded domains. As analytical solutions are available only for simple geometries and material properties, their application is limited [4].The approximate methods are, in general, spatially and temporally local. They are computationally efficient by themselves, but have to be applied at a so-called artificial boundary sufficiently far away from the region of interest in order to obtain results of acceptable accuracy. This increases the total computational effort. The viscous boundary [15] consists of dashpots that absorb plane waves propagating perpendicularly to the artificial boundary. The Bayliss, Gunzburger and Turkel (BGT) boundary [16] annihilates selected terms of cylindrical or spherical waves. Infinite elements [5, 17] use decay functions representing the wave propagation toward infinity as shape functions of the displacements. In the method of perfectly matched layer [18], the domain of interest is surrounded by a fini...