An exact formula is obtained for the term of order c2 in the expression for the bulk stress in a suspension of force-free spherical particles in Newtonian ambient fluid, where c is the volume fraction of the spheres and c [Lt ] 1. The particles may be of different sizes, and composed of either solid or fluid of arbitrary viscosity. The method of derivation circumvents the familiar obstacle, of non-absolutely convergent integrals representing the effect of all pair interactions in which one specified particle takes part, by the judicious use of a certain quantity which is affected by the presence of distant particles in a similar way and whose mean value is known exactly. The bulk stress is in general of non-Newtonian form and depends on the statistical properties of the suspension which in turn are dependent on the type of bulk flow.The formula contains two functions which are parameters of the flow field due to two spherical particles immersed in fluid in which the velocity gradient is uniform at infinity. One of them, p(r, t), represents the probability density for the vector r separating the centres of the two particles. The variation of p(r, t) for a moving material point in r-space due to hydrodynamic action is found in terms of a function q(r), and this gives p(r, t) explicitly over the whole of the region of r-space occupied by trajectories of one particle centre relative to another which come from infinity. In a region of closed trajectories, steady-state hydrodynamic action alone does not determine the relation between the values of p (r, t) for different material points. The function q(r) is singular when the spheres touch, and the contribution of nearly-touching spheres to the bulk stress is evidently important. Approximate numerical values of all the relevant functions are presented for the case of rigid spherical particles of uniform size.In the case of steady pure straining motion of the suspension, all trajectories in r-space come from infinity, the suspension has isotropic structure and the stress behaviour can be represented (to order c2) in terms of an effective viscosity . It is estimated from the available numerical data that for a suspension of identical rigid spherical particles \[ {\mathop\mu\limits^{*}}/\mu = 1 + 2.5c + 7.6c^2, \] the error bounds on the coefficient of c2 being about ∓ 0.8. In the important case of steady simple shearing motion, there is a region of closed trajectories of one sphere centre relative to another, of infinite volume. The stress system is here not of Newtonian form, and numerical results are not obtainable until the probability, density p(r, t) can be made determinate in the region of closed trajectories by the introduction of some additional physical process, such as three-sphere encounters or Brownian motion, or by the assumption of some particular initial state.In the analogous problem for an incompressible solid suspension it may be appropriate to assume that for many methods of manufacture p(r, t) is uniform over the accessible part of r-space, in which event the solid suspension has ‘Newtonian’ elastic behaviour and the ratio of the effective shear modulus to that of the matrix is estimated to be 1 + 2·5c + 5·2c2 for a suspension of identical rigid spheres.