The finite size of the simulated system may have a significant impact on the location of phase boundaries determined from a simulation with periodic boundary conditions. Only in the limit of infinite system size, such simulations converge to the properties of bulk phase-separated systems. To demonstrate that our results were not significantly affected by the size of the simulated system, in Sec. S2.3 in the ESI we present structural analysis of systems with various sizes, determined by the numbers of polymer chains in the polyelectrolyte complex, N A = N C , ranging from 32 up to 96. In a typical simulation, we observe that L ≈ 8R g , where R g is the mean radius of gyration of the polyions. Therefore the box is large enough, so that chains do not interact strongly with their own periodic images. Accordingly, we observe in Sec. S2.3 that the positions of binodals and partition coefficients of the ions are virtually the same for all simulated sizes of the system. Therefore, we used N A = N C = 32 in all further simulations, which had the most favourable computational costs.