The study of axially symmetric stationary multi-black-hole configurations and the force between co-axially rotating black holes involves, as a first step, an analysis on the “boundary regularity” of the so-called reduced singular harmonic maps. We carry out this analysis by considering those harmonic maps as solutions to some homogeneous divergence systems of partial differential equations with singular coefficients. Our results extend previous works by Weinstein (Comm Pure Appl Math 43:903–948, 1990; Comm Pure Appl Math 45:1183–1203, 1992) and by Li and Tian (Manu Math 73(1):83–89, 1991; Commun Math Phys 149:1–30, 1992; Differential geometry: PDE on manifolds, vol 54, pp. 317–326, 1993). This paper is based on the Ph.D. thesis of the author (Singular harmonic maps into hyperbolic spaces and applications to general relativity, PhD thesis, The State University of New Jersey, Rutgers, 2009).