We investigate the strongly correlated ion dynamics and the degree of coupling achievable in the evolution of freely expanding ultracold neutral plasmas. We demonstrate that the ionic Coulomb coupling parameter Γi increases considerably in later stages of the expansion, reaching the strongly coupled regime despite the well-known initial drop of Γi to order unity due to disorder-induced heating. Furthermore, we formulate a suitable measure of correlation and show that Γi calculated from the ionic temperature and density reflects the degree of order in the system if it is sufficiently close to a quasisteady state. At later times, however, the expansion of the plasma cloud becomes faster than the relaxation of correlations, and the system does not reach thermodynamic equilibrium anymore.PACS numbers: 52.27. Gr,32.80.Pj,05.70.Ln Freely expanding ultracold neutral plasmas (UNPs) [1] have attracted wide attention both experimentally [2,3,4,5] and theoretically [6,7,8,9,10]. A main motivation of the early experiments was the creation of a strongly coupled plasma, with the Coulomb coupling parameter (CCP) Γ = e 2 /(ak B T ) ≫ 1 (where T is temperature and a is the Wigner-Seitz radius). From the experimental setup of [1], the CCPs of electrons and ions were estimated to be of the orders of Γ e ≈ 30 and Γ i ≈ 30000, respectively. By changing the frequency of the ionizing laser, the electronic temperature can be varied, offering the prospect of controlling the coupling strength of the electrons and creating UNPs where either one, namely the ionic, or both components could be strongly coupled.However, due to unavoidable heating effects [6,11,12] these hopes have not materialized yet, and only Γ e ≈ 0.2 and Γ i ≈ 2 have been confirmed. Furthermore, the evolution of the expanding plasma turns out to be a rather intricate problem of non-equilibrium plasma physics for which a clear definition of the degree of correlation is not obvious to begin with.The goal of this letter is twofold: Firstly, we will formulate a consistent measure of correlation for expanding ultracold plasmas, and secondly we demonstrate that the strongly correlated regime with Γ i ≈ 10 for the ionic plasma component can be reached by simply waiting until the plasma has (adiabatically) expanded long enough under already realized experimental conditions. This is remarkable in the light of alternatives proposed to increase Γ i [12,13,14,15,16,17] which are experimentally rather involved.Substantiating both of our statements theoretically requires the ability to propagate the plasma numerically over a long time with full account of the ionic correlations. To this end, we have developed a hybrid molecular dynamics (H-MD) method [9] for the description of ultracold neutral plasmas. In our approach, ions and recombined atoms are propagated in the electronic mean-field potential with the full ion-ion interaction taken into account. The much faster and weakly coupled electrons, on the other hand, are treated on a hydrodynamical level. Elastic as well as inelastic ...