This white paper describes the science case for Very Long Baseline Interferometry (VLBI) and provides suggestions towards upgrade paths for the European VLBI Network (EVN). The EVN is a distributed long-baseline radio interferometric array, that operates at the very forefront of astronomical research. Recent results, together with the new science possibilities outlined in this vision document, demonstrate the EVN's potential to generate new and exciting results that will transform our view of the cosmos. Together with e-MERLIN, the EVN provides a range of baseline lengths that permit unique studies of faint radio sources to be made over a wide range of spatial scales. The science cases are reviewed in six chapters that cover the following broad areas: cosmology, galaxy formation and evolution, innermost regions of active galactic nuclei, explosive phenomena and transients, stars and stellar masers in the Milky Way, celestial reference frames and space applications. The document concludes with identifying the synergies with other radio, as well as multi-band/multi-messenger instruments, and provide the recommendations for future improvements. The appendices briefly describe other radio VLBI arrays, the technological framework for EVN developments, and a selection of spectral lines of astrophysical interest below 100 GHz. The document includes a glossary for non-specialists, and a list of acronyms at the end.

We present a study of the stellar host galaxy, CO (1$-$0) molecular gas distribution and AGN emission on 50 to 500 pc-scales of the gravitationally lensed dust-obscured AGN MG J0751+2716 and JVAS B1938+666 at redshifts 3.200 and 2.059, respectively. By correcting for the lensing distortion using a grid-based lens modelling technique, we spatially locate the different emitting regions in the source plane for the first time. Both AGN host galaxies have 300 to 500 pc-scale size and surface brightness consistent with a bulge/pseudo-bulge, and 2 kpc-scale AGN radio jets that are embedded in extended molecular gas reservoirs that are 5 to 20 kpc in size. The CO (1$-$0) velocity fields show structures possibly associated with discs (elongated velocity gradients) and interacting objects (off-axis velocity components). There is evidence for a decrement in the CO (1$-$0) surface brightness at the location of the host galaxy, which may indicate radiative feedback from the AGN, or offset star formation.We find CO-H$_2$ conversion factors of around $\alpha_{\rm CO} = 1.5\pm0.5$ (K km s$^{-1}$ pc$^2$)$^{-1}$, molecular gas masses of $> 3\times10^{10}$ M$_{\odot}$, dynamical masses of $\sim 10^{11}$ M$_{\odot}$ and gas fractions of around 60 per cent. The intrinsic CO line luminosities are comparable to those of unobscured AGN and dusty star-forming galaxies at similar redshifts, but the infrared luminosities are lower, suggesting that the targets are less efficient at forming stars. Therefore, they may belong to the AGN feedback phase predicted by galaxy formation models, because they are not efficiently forming stars considering their large amount of molecular gas.

We present a new gravitational lens modelling technique designed to model high-resolution interferometric observations with large numbers of visibilities without the need to pre-average the data in time or frequency. We demonstrate the accuracy of the method using validation tests on mock observations. Using small data sets with $\sim 10^3$ visibilities, we first compare our approach with the more traditional direct Fourier transform (DFT) implementation and direct linear solver. Our tests indicate that our source inversion is indistinguishable from that of the DFT. Our method also infers lens parameters to within 1 to 2 per cent of both the ground truth and DFT, given sufficiently high signal-to-noise ratio (SNR). When the SNR is as low as 5, both approaches lead to errors of several tens of per cent in the lens parameters and a severely disrupted source structure, indicating that this is an issue related to the data quality rather than the modelling technique of choice. We then analyze a large data set with $\sim 10^8$ visibilities and a SNR matching real global Very Long Baseline Interferometry observations of the gravitational lens system MG J0751+2716. The size of the data is such that it cannot be modelled with traditional implementations. Using our novel technique, we find that we can infer the lens parameters and the source brightness distribution, respectively, with an RMS error of 0.25 and 0.97 per cent relative to the ground truth.

We present a multi-wavelength analysis of two highly magnified strong gravitationally lensed galaxies, CLASS B0712+472 and CLASS B1608+656, at redshifts $1.34$ and $1.394$, respectively, using new VLBI and archival HST observations. We reconstruct the positions of the radio and optical emissions with their uncertainties using Monte Carlo sampling. We find that in CLASS B0712+472 the optical and radio emissions are co-spatial within $2\pm5$ mas ($17\pm 42$ pc at redshift of 1.34). But, in CLASS B1608+656, we reconstruct an optical-radio offset of $25\pm16$ mas ($214\pm137$ pc at redshift of 1.394), the smallest offset measured for an AGN at such high redshift. The spectral features indicate that CLASS B1608+656 is a post-merger galaxy, which, in combination with the optical-VLBI offset reported here, makes CLASS B1608+656 a promising candidate for a high-$z$ offset-AGN. Furthermore, the milliarcsecond angular resolution of the VLBI observations combined with the precise lens models allow us to spatially locate the radio emission at $0.05$ mas precision ($0.4$ pc) in CLASS B0712+472, and $0.009$ mas precision ($0.08$ pc) in CLASS B1608+656. The search for optical-radio offsets in high redshift galaxies will be eased by the upcoming synoptic all-sky surveys, including E-ELT and SKA, which are expected to find $\sim 10^5$ strongly lensed galaxies, opening an era of large strong lensing samples observed at high angular resolution.

We present the discovery of PSO J030947.49+271757.31, the radio brightest (23.7 mJy at 1.4 GHz) active galactic nucleus (AGN) at z>6.0. It was selected by cross-matching the NRAO VLA Sky Survey and the Panoramic Survey Telescope and Rapid Response System PS1 databases and its high-z nature was confirmed by a dedicated spectroscopic observation at the Large Binocular Telescope. A pointed Neil Gehrels $Swift$ Observatory XRT observation allowed us to measure a flux of $\sim$3.4$\times$10$^{-14}$ erg s$^{-1}$ cm$^{-2}$ in the [0.5-10] keV energy band, which also makes this object the X-ray brightest AGN ever observed at z>6.0. Its flat radio spectrum ($\alpha_{\nu r}$<0.5), very high radio loudness (R>10$^3$), and strong X-ray emission, compared to the optical, support the hypothesis of the blazar nature of this source. %i.e. a radio-loud (RL) AGN with the relativistic jet pointed toward us. Assuming that this is the only blazar at this redshift in the surveyed area of sky, we derive a space density of blazars at z$\sim$6 and with M$_{1450 \mbox{\AA}}$ < -25.1 of 5.5$^{+11.2}_{-4.6}$$\times$10$^{-3}$ Gpc$^{-3}$. From this number, and assuming a reasonable value of the bulk velocity of the jet ($\Gamma$=10), we can also infer a space density of the entire radio-loud AGN population at z$\sim$6 with the same optical/UV absolute magnitude of 1.10$^{+2.53}_{-0.91}$ Gpc$^{-3}$. Larger samples of blazars will be necessary to better constrain these estimates.

In this paper, we exploit the gravitational lensing effect to detect proper motion in the highly magnified gravitationally lensed source MG B2016+112. We find positional shifts up to 6 mas in the lensed images by comparing two Very Long Baseline Interferometric (VLBI) radio observations at 1.7 GHz that are separated by 14.359 years, and provide an astrometric accuracy of the order of tens of $\mu$as. From lens modelling, we exclude a shift in the lensing galaxy as the cause of the positional change of the lensed images, and we assign it to the background source. The source consists of four sub-components separated by $\sim 175$ pc, with proper motion of the order of tens $\mu$as yr$^{-1}$ for the two components at highest magnification ($\mu\sim350$) and of the order of a few mas yr$^{-1}$ for the two components at lower magnification ($\mu\sim2$). We propose single AGN and dual AGN scenarios to explain the source plane. Although, the latter interpretation is supported by the archival multi-wavelength properties of the object. In this case, MG B2016+112 would represent the highest redshift dual radio-loud AGN system discovered thus far, and would support the merger interpretation for such systems. Also, given the low probability ($\sim10^{-5}$) of detecting a dual AGN system that is also gravitationally lensed, if confirmed, this would suggest that such dual AGN systems must be more abundant in the early Universe than currently thought.

Gravitational lensing is a powerful tool for quantifying the mass content and distribution in distant galaxies. By using milliarcsecond angular resolution observations of radio-loud gravitationally lensed sources it is also possible to detect and quantify small deviations from a smooth mass density distribution, which can be due to low mass substructures in the lensing galaxy. We present high-resolution global VLBI observations of the gravitationally lensed radio source MG J0751+2716 (at z = 3.2), that shows evidence of both compact and extended structure (core-jet morphology) across several gravitational arcs. These data provide a wealth of observational constraints that are used to determine the inner (baryonic and dark matter) mass profile of a group of galaxies and also investigate the smoothness of the dark matter distribution on mas-scales, which is sensitive to possible structures of $10^{6-7}$ M$_{\odot}$ within the lensing halo or along the line-of-sight. Our lens modelling finds evidence for astrometric anomalies in this system, which suggest presence of extra mass structure in the lens model. To date this kind of detailed studies of gravitational lensing systems like MG J0751+2716 has been limited by the currently small sample of radio-loud gravitational lenses. In this context, we also present a new pilot gravitational lens search in the VLBI survey mJIVE-20, in perspective of future surveys with the next generation of radio interferometers.

We present a novel pilot search for gravitational lenses in the mJIVE-20 survey, which observed $24\,903$ radio sources selected from FIRST with the VLBA at an angular resolution of 5 mas. We have taken the visibility data for an initial $3\,640$ sources that were detected by the mJIVE-20 observations and re-mapped them to make wide-field images, selecting fourteen sources that had multiple components separated by $\geq100$ mas, with a flux-ratio of $\leq15$:$1$ and a surface brightness consistent with gravitational lensing. Two of these candidates are re-discoveries of gravitational lenses found as part of CLASS. The remaining twelve candidates were then re-observed at 1.4 GHz and then simultaneously at 4.1 and 7.1 GHz with the VLBA to measure the spectral index and surface brightness of the individual components as a function of frequency. Ten were rejected as core-jet or core-hotspot(s) systems, with surface brightness distributions and/or spectral indices inconsistent with gravitational lensing, and one was rejected after lens modelling demonstrated that the candidate lensed images failed the parity test. The final lens candidate has an image configuration that is consistent with a simple lens mass model, although further observations are required to confirm the lensing nature. Given the two confirmed gravitational lenses in the mJIVE-20 sample, we find a robust lensing-rate of $1$:($318\pm225$) for a statistical sample of 635 radio sources detected on mas-scales, which is consistent with that found for CLASS.

We present milliarcsecond (mas) angular resolution observations of the gravitationally lensed radio source MG J0751+2716 (at z=3.2) obtained with global Very Long Baseline Interferometry (VLBI) at 1.65 GHz. The background object is highly resolved in the tangential and radial directions, showing evidence of both compact and extended structure across several gravitational arcs that are 200 to 600~mas in size. By identifying compact sub-components in the multiple images, we constrain the mass distribution of the foreground z=0.35 gravitational lens using analytic models for the main deflector [power-law elliptical mass model; $\rho(r) \propto r^{-\gamma}$, where $\gamma=2$ corresponds to isothermal] and for the members of the galaxy group. Moreover, our mass models with and without the group find an inner mass-density slope steeper than isothermal for the main lensing galaxy, with $\gamma_1 = 2.08 \pm 0.02$ and $\gamma_2 = 2.16 \pm 0.02$ at the 4.2$\sigma$ level and 6.8$\sigma$ level, respectively, at the Einstein radius ($b_1 = 0.4025 \pm 0.0008$ and $b_2 = 0.307 \pm 0.002$ arcsec, respectively). We find randomly distributed image position residuals of about 3 mas, which are much larger that the measurement errors ($40$ $\mu$as on average). This suggests that at the mas level, the assumption of a smooth mass distribution fails, requiring additional structure in the model. However, given the environment of the lensing galaxy, it is not clear whether this extra mass is in the form of sub-haloes within the lens or along the line of sight, or from a more complex halo for the galaxy group.

Flux ratio anomalies in quasar lenses can be attributed to dark matter substructure surrounding the lensing galaxy and, thus, used to constrain the substructure mass fraction. Previous applications of this approach infer a substructure abundance that potentially in tension with the predictions of a $\Lambda$CDM cosmology. However, the assumption that all flux ratio anomalies are due to substructure is a strong one, and alternative explanations have not been fully investigated. Here, we use new high-resolution near-IR Keck~II adaptive optics imaging for the lens system CLASS B0712+472 to perform pixel-based lens modelling for this system and, in combination with new VLBA radio observations, show that the inclusion of the disc in the lens model can explain the flux ratio anomalies without the need for dark matter substructures. The projected disc mass comprises 16% of the total lensing mass within the Einstein radius and the total disc mass is $1.79 \times 10^{10} M_{sun}$. The case of B0712+472 adds to the evidence that not all flux ratio anomalies are due to dark subhaloes, and highlights the importance of taking the effects of baryonic structures more fully into account in order to obtain an accurate measure of the substructure mass fraction.

Gravitational lens flux-ratio anomalies provide a powerful technique for measuring dark matter substructure in distant galaxies. However, before using these flux-ratio anomalies to test galaxy formation models, it is imperative to ascertain that the given anomalies are indeed due to the presence of dark matter substructure and not due to some other component of the lensing galaxy halo or to propagation effects. Here we present the case of CLASS~B1555+375, which has a strong radio-wavelength flux-ratio anomaly. Our high-resolution near-infrared Keck~II adaptive optics imaging and archival Hubble Space Telescope data reveal the lensing galaxy in this system to have a clear edge-on disc component that crosses directly over the pair of images that exhibit the flux-ratio anomaly. We find that simple models that include the disc can reproduce the cm-wavelength flux-ratio anomaly without requiring additional dark matter substructure. Although further studies are required, our results suggest the assumption that all flux-ratio anomalies are due to a population of dark matter sub-haloes may be incorrect, and analyses that do not account for the full complexity of the lens macro-model may overestimate the substructure mass fraction in massive lensing galaxies.

We present results on multifrequency Very Long Baseline Array (VLBA) monitoring observations of the double-image gravitationally lensed blazar JVAS B0218+357. Multi-epoch observations started less than one month after the gamma-ray flare detected in 2012 by the Large Area Telescope on board Fermi, and spanned a 2-month interval. The radio light curves did not reveal any significant flux density variability, suggesting that no clear correlation between the high energy and low-energy emission is present. This behaviour was confirmed also by the long-term Owens Valley Radio Observatory monitoring data at 15 GHz. The milliarcsecond-scale resolution provided by the VLBA observations allowed us to resolve the two images of the lensed blazar, which have a core-jet structure. No significant morphological variation is found by the analysis of the multi-epoch data, suggesting that the region responsible for the gamma-ray variability is located in the core of the AGN, which is opaque up to the highest observing frequency of 22 GHz.