The Gross-Pitaevskii equation -which describes interacting bosons in the mean-field approximation -possesses solitonic solutions in dimension one. For repulsively interacting particles, the stationary soliton is dark, i.e. is represented by a local density minimum. Many-body effects may lead to filling of the dark soliton. Using quasi-exact many-body simulations, we show that, in single realizations, the soliton appears totally dark although the single particle density tends to be uniform. [4][5][6][7][8][9]. Within the mean field approach and at zero temperature, an atomic Bose-Einstein condensate (BEC) is described by the Gross-Pitaevskii equation (GPE). In one-dimensional (1D) space, this non-linear wave equation possesses a dark (bright) soliton solution if particle interactions are effectively repulsive (attractive). The GPE assumes all particles in the same single particle state and neglects that the interactions can populate other modes.There is a debate in the literature concerning manybody effects in a dark soliton [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. The GPE predicts a stable dark soliton state in 1D. The Bogoliubov corrections show that signatures of the soliton tend to disappear in the single particle density because particles depleted from the condensate fill the soliton notch [12][13][14][15]. However, the single particle density refers to results averaged over many experimental realizations and is not able to predict a single experimental outcome. Within the Bogoliubov approach, it is possible to simulate single experiments: in a single photo, the soliton is completely dark but its position varies randomly from realization to realization [14,15]. This is a signature of the quantum character of the soliton, i.e. the soliton position is not described by a classical variable but by a probability density. The Bogoliubov analysis is limited to quantum fluctuations of the soliton position on a scale smaller than the healing length and cannot predict what happens when many-body effects become stronger. The latter situation has been considered in Refs. [19,20] within many-body numerical calculations. Analysis of the second order correlation function leads the authors to the conclusion that a dark soliton cannot be observed when the single particle density becomes uniform. However, full simulations of single experiments involve much higher order correlation functions [21,22]. It is the goal of the present work to perform many-body numerical simulations of the entire experiment starting from soliton preparation till the destructive measurement of all particle positions. We show that a fully dark soliton can be observed when manybody effects are strong, although its measured position is random and varies among experimental realizations. In the following we assume the zero temperature limit and do not consider thermodynamical instability of the dark soliton [26,27].The stationary solution of the 1D GPE for particles of mass m and interaction coefficient g 0 , associated with a da...