Gaia kinematics reveal a complex lopsided and twisted Galactic disc warp

There are few warp kinematic models able to characterise both structure and kinematics. These models are needed to shed some light to the lopsidedness and twisting of the stellar warp already seen in gas and dust. We use the vertical information from the Gaia Data Release 2 up to G=20 mag to study the structure of the Galactic warp, the related vertical motions and the relation with age. We analyse two populations up to galactocentric distances of R=16 kpc, a young bright sample formed by OB stars and an older one of Red Giant Branch stars. We use two methods (ngc PCM and LonKin) based on the Gaia observables, together with 2D maps of positions and proper motions in the Galactic plane. We indeed confirm the age dependency of the Galactic warp, both in positions and kinematics, being the height of the Galactic warp ~ 0.2 kpc for the OB sample and ~1.0 kpc for the Red Giant Branch at R=14 kpc. Both methods find that the onset radius is 12 10 kpc the line-of-nodes twists away from the Sun-anticentre line towards galactic azimuths ~180-200 deg increasingly with radius. Also, the RGB sample reveals a slightly lopsided stellar warp with |Z(down)|- |Z(up)|~250 pc. The line of maximum of proper motions in latitude is systematically offset from the line-of-nodes estimated from the spatial data, which our warp models predict as a kinematic signature of lopsidedness. We also show a prominent wave-like pattern of a bending mode different in the OB and RGB samples. Both positions and kinematics also reveal substructures that might not be related to the large scale Galactic warp or to the bending mode. GaiaDR2 reveals a high degree of complexity that triggers the need for complex kinematic models, flexible enough to combine wave-like patterns and an S-shaped lopsided and twisted warp [abridged].