Data from a dedicated cosmic ray run of the ALEPH detector were used in a study of muon trident production, i.e., muon pairs produced by muons. Here the overburden and the calorimeters are the target materials while the ALEPH time projection chamber provides the momentum measurements. A theoretical estimate of the muon trident cross section is obtained by developing a Monte Carlo simulation for muon propagation in the overburden and the detector. Two muon trident candidates were found to match the expected theoretical pattern. The observed production rate implies that the nuclear form factor cannot be neglected for muon tridents. DOI: 10.1103/PhysRevLett.96.021801 PACS numbers: 13.40.Gp, 14.60.Ef, 25.30.Mr, 95.85.Ry The trident interaction process occurs when a lepton l 1 produces a pair of leptons l 2 ; l ÿ 2 in the field of nuclear charge Ze, thus l 1 Z; l 2 l ÿ 2 Zl 1 . Good theoretical estimates of the trident production rate date back to 1937 [1,2]. The first papers on this subject considered: small scattering angles, no exchange between diagrams [1], and only part of the leading order Feynman diagrams. The channels e e ÿ and ÿ were at this level indistinguishable except for the mass difference. About 30 years later the full QED formula of the differential trident cross section was computed [3][4][5]. However, there were two main difficulties to be solved: inclusion of the atomic and nuclear form factors and the simultaneous evaluation of the phase-space integral. Because of the high complexity of the full integral formula, the obtained QED result is impractical for most cross-section estimations. In practical applications, the result of a numerical integration of the QED matrix element over a multidimensional phase space is approximated by an analytical parametrization which is convenient for Monte Carlo (MC) simulations of leptonic tridents.The early studies were motivated: by the necessity to compute radiative energy losses for charged particles (e.g., muon energy loss by emission of e e ÿ in the field of a nucleus), by the unknown nature of the muon [fermion or boson, see [6], concluded in [7] ], and by possible anomalous muon-muon interactions for Z; ÿ Z [8]. Interference between amplitudes of different Feynman diagrams was estimated to be negligible for muons with energies higher than 10 GeV [9]. Because the study presented in this Letter concerns high energy muons and due to the decreasing probability of muon tridents at low energy, interference terms were ignored.In 1967 Kelner [10] computed the QED differential cross section for tridents neglecting the nuclear form factor. The expression accounts for both Z; ÿ Z and Z; e e ÿ Z and its subsequent parametrization is the Bugaev-Kotov-Rozental formula [11] [BKR in [12] ] that reduces the multidimensional integral to an analytical result. The formula has become a standard reference for professional muon transport Monte Carlo programs [e.g., [13] ] and for standard muon energy loss tables [14]. The BKR formula was later corrected by Kokoulin, Kelner ...