We analyze the time-dependent formation of the spectral function of an Anderson impurity model in the Kondo regime within a numerically exact real-time quantum Monte Carlo framework. At steady state, splitting of the Kondo peak occurs with nontrivial dependence on voltage and temperature, and with little effect on the location or intensity of high-energy features. Examining the transient development of the Kondo peak after a quench from an initially uncorrelated state reveals a two-stage process where the initial formation of a single central Kondo peak is followed by splitting. We analyze the time dependence of splitting in detail and demonstrate a strong dependence of its characteristic timescale on the voltage. We expect both the steady state and the transient phenomenon to be experimentally observable.Introduction. Interacting quantum many-body systems often exhibit highly entangled states that cannot be described within an independent particle formalism. The Kondo effect in a quantum dot 1,2 coupled to noninteracting leads is the paradigmatic example for such a state, as the dot electrons hybridize with the leads to form a highly correlated Kondo singlet state 3 . This state manifests itself as a sharp peak in the local density of states 2,4 . The establishment of Kondo correlations can be examined in a quantum quench scenario, where an initially uncorrelated state slowly develops a coherence peak over time 5,6 .In the presence of a voltage, the Kondo peak is strongly suppressed and splits into two smaller peaks 7-11 . Previous work has argued that the peak-to-peak distance is given by the voltage 12-17 and that the split state is significantly less correlated than the equilibrium state 12 . It is therefore natural to examine the establishment of splitting after a quench from an initially uncorrelated state, and to expect that this less correlated state forms on a timescale shorter than that of the equilibrium state.Despite significant analytical progress 18-24 , an accurate investigation of this scenario requires numerical methods that are able to simulate the real-time evolution after a quench accurately, for times long enough to reach the steady state. Additionally, a full account of the continuous lead spectrum is crucial for correct treatment of the nonequilibrium steady state. The major families of numerical methods include the noncrossing approximation and its higher-order generalizations 25 , wavefunction-based methods 26-31 , real-time path integral techniques 32-35 , the time-dependent numerical renormalization group 36-40 , hierarchical equations of motion 41-44 , the auxiliary master equation approach 45-49 , and a wide variety of quantum Monte Carlo methods 50-66 . Most of these approaches fall short in at least one of the aforementioned requirements. This situation has changed with the development of the numerically exact inchworm quantum Monte Carlo method 67-71 that in many cases eliminates the dynamical sign problem and is thereby able to reach the relevant timescales.