The Central Molecular Zone (CMZ) is the most extreme star-forming environment in the Milky Way in terms of gas pressures, densities, temperatures, and dynamics. It acts as a critical test bed for developing star formation theories applicable to the (high-redshift-like) conditions under which most stars in the Universe formed. We present a set of numerical simulations of molecular clouds orbiting on the 100-pc stream that dominates the molecular gas reservoir of the CMZ, with the goal of characterising their morphological and kinematic evolution in response to the external gravitational potential and their eccentric orbital motion. These simulations capture the evolution of single clouds in a strong and plausibly dominant background potential. We find that the evolution of the clouds is closely coupled to the orbital dynamics and their arrival on the 100-pc stream marks a transformative event in their lifecycle. The clouds' sizes, aspect ratios, position angles, filamentary structure, column densities, velocity dispersions, line-of-sight velocity gradients, angular momenta, and overall kinematic complexity are controlled by the background potential and their passage through the orbit's pericentre. We compare these predictions of our simulations to observations of clouds on the Galactic Centre `dust ridge' and find that the inclusion of galactic dynamics naturally reproduces a surprisingly wide variety of key observed morphological and kinematic features. We argue that the accretion of gas clouds onto the central regions of galaxies, where the rotation curve turns over and the tidal field becomes fully compressive, is likely to lead to their collapse and associated star formation. This can generate an evolutionary progression of cloud collapse with a common zero point. Together, these processes may naturally give rise to the synchronised starbursts observed in numerous galactic nuclei.

We present cosmological hydrodynamical simulations including atomic and molecular non-equilibrium chemistry, multi-frequency radiative transfer (0.7-100 eV sampled over 150 frequency bins) and stellar population evolution to investigate the host candidates of the seeds of supermassive black holes coming from direct collapse of gas in primordial haloes (direct-collapse black holes, DCBHs). We consistently address the role played by atomic and molecular cooling, stellar radiation and metal spreading of C, N, O, Ne, Mg, Si, S, Ca, Fe, etc. from primordial sources, as well as their implications for nearby quiescent proto-galaxies under different assumptions for early source emissivity, initial mass function and metal yields. We find that putative DCBH host candidates need powerful primordial stellar generations, since common solar-like stars and hot OB-type stars are neither able to determine the conditions for direct collapse nor capable of building up a dissociating Lyman-Werner background radiation field. Thermal and molecular features of the identified DCBH host candidates in the scenario with very massive primordial stars seem favourable, with illuminating Lyman-Werner intensities featuring values of 1-50 J21. Nevertheless, additional non-linear processes, such as merger events, substructure formation, rotational motions and photo-evaporation, should inhibit pure DCBH formation in 2/3 of the cases. Local turbulence may delay gas direct collapse almost irrespectively from other environmental conditions. The impact of large Lyman-Werner fluxes at distances smaller than 5 kpc is severely limited by metal pollution.

The evolution of molecular clouds in galactic centres is thought to differ from that in galactic discs due to a significant influence of the external gravitational potential. We present a set of numerical simulations of molecular clouds orbiting on the 100-pc stream of the Central Molecular Zone (the central $\sim500$ pc of the Galaxy) and characterise their morphological and kinematic evolution in response to the background potential and eccentric orbital motion. We find that the clouds are shaped by strong shear and torques, by tidal and geometric deformation, and by their passage through the orbital pericentre. Within our simulations, these mechanisms control cloud sizes, aspect ratios, position angles, filamentary structure, column densities, velocity dispersions, line-of-sight velocity gradients, spin angular momenta, and kinematic complexity. By comparing these predictions to observations of clouds on the Galactic Centre 'dust ridge', we find that our simulations naturally reproduce a broad range of key observed morphological and kinematic features, which can be explained in terms of well-understood physical mechanisms. We argue that the accretion of gas clouds onto the central regions of galaxies, where the rotation curve turns over and the tidal field is fully compressive, is accompanied by transformative dynamical changes to the clouds, leading to collapse and star formation. This can generate an evolutionary progression of cloud collapse with a common starting point, which either marks the time of accretion onto the tidally-compressive region or of the most recent pericentre passage. Together, these processes may naturally produce the synchronised starbursts observed in numerous (extra)galactic nuclei.

We investigate four young, but non-accreting, very low mass stars in Orion, which show irregular eclipses by circumstellar dust. The eclipses are not recurring periodically, are variable in depth, lack a flat bottom, and their duration is comparable to the typical timescale between eclipses. The dimming is associated with reddening consistent with dust extinction. Taken together this implies the presence of rings around these four stars, with radii ranging from 0.01 to 40 AU, comprised of optically thin dust clouds. The stars also show IR excess indicating the presence of evolved circumstellar disks, with orders of magnitude more material than needed for the eclipses. However, the rings need to cover an opening angle of about 20 degrees to explain how common these variable stars are in the coeval population in the same region, which is more extended than a typical disk. Thus, we propose that the rings may not be part of the disks, but instead separate structures with larger scale heights. To be sustained over years, the rings need to be replenished by dust from the disk or gravitationally bound to an object (e.g., planets or planetesimals). These four stars belong to a growing and diverse class of post-T Tauri stars with dips or eclipses in their lightcurves. Dusty rings with scale heights exceeding those of disks may be a common phenomenon at stellar ages between 5 and 10 Myr, in the transition from accretion disks to debris disks. These structures could be caused by migrating planets and may be signposts for the presence of young planetary systems.

This article describes a data center hosting a web portal for accessing and sharing the output of large, cosmological, hydro-dynamical simulations with a broad scientific community. It also allows users to receive related scientific data products by directly processing the raw simulation data on a remote computing cluster. The data center has a multi-layer structure: a web portal, a job control layer, a computing cluster and a HPC storage system. The outer layer enables users to choose an object from the simulations. Objects can be selected by visually inspecting 2D maps of the simulation data, by performing highly compounded and elaborated queries or graphically by plotting arbitrary combinations of properties. The user can run analysis tools on a chosen object. These services allow users to run analysis tools on the raw simulation data. The job control layer is responsible for handling and performing the analysis jobs, which are executed on a computing cluster. The innermost layer is formed by a HPC storage system which hosts the large, raw simulation data. The following services are available for the users: (I) {\sc ClusterInspect} visualizes properties of member galaxies of a selected galaxy cluster; (II) {\sc SimCut} returns the raw data of a sub-volume around a selected object from a simulation, containing all the original, hydro-dynamical quantities; (III) {\sc Smac} creates idealised 2D maps of various, physical quantities and observables of a selected object; (IV) {\sc Phox} generates virtual X-ray observations with specifications of various current and upcoming instruments.

Voronoi grids have been successfully used to represent density structures of gas in astronomical hydrodynamics simulations. While some codes are explicitly built around using a Voronoi grid, others, such as Smoothed Particle Hydrodynamics (SPH), use particle-based representations and can benefit from constructing a Voronoi grid for post-processing their output. So far, calculating the density of each Voronoi cell from SPH data has been done numerically, which is both slow and potentially inaccurate. This paper proposes an alternative analytic method, which is fast and accurate. We derive an expression for the integral of a cubic spline kernel over the volume of a Voronoi cell and link it to the density of the cell. Mass conservation is ensured rigorously by the procedure. The method can be applied more broadly to integrate a spherically symmetric polynomial function over the volume of a random polyhedron.

We present an implementation of thermal conduction including the anisotropic effects of magnetic fields for SPH. The anisotropic thermal conduction is mainly proceeding parallel to magnetic fields and suppressed perpendicular to the fields. We derive the SPH formalism for the anisotropic heat transport and solve the corresponding equation with an implicit conjugate gradient scheme. We discuss several issues of unphysical heat transport in the cases of extreme ansiotropies or unmagnetized regions and present possible numerical workarounds. We implement our algorithm into the GADGET code and study its behaviour in several test cases. In general, we reproduce the analytical solutions of our idealised test problems, and obtain good results in cosmological simulations of galaxy cluster formations. Within galaxy clusters, the anisotropic conduction produces a net heat transport similar to an isotropic Spitzer conduction model with an efficiency of one per cent. In contrast to isotropic conduction our new formalism allows small-scale structure in the temperature distribution to remain stable, because of their decoupling caused by magnetic field lines. Compared to observations, isotropic conduction with more than 10 per cent of the Spitzer value leads to an oversmoothed temperature distribution within clusters, while the results obtained with anisotropic thermal conduction reproduce the observed temperature fluctuations well. A proper treatment of heat transport is crucial especially in the outskirts of clusters and also in high density regions. It's connection to the local dynamical state of the cluster also might contribute to the observed bimodal distribution of cool core and non cool core clusters. Our new scheme significantly advances the modelling of thermal conduction in numerical simulations and overall gives better results compared to observations.

In spring 2015, the Leibniz Supercomputing Centre (Leibniz-Rechenzentrum, LRZ), installed their new Peta-Scale System SuperMUC Phase2. Selected users were invited for a 28 day extreme scale-out block operation during which they were allowed to use the full system for their applications. The following projects participated in the extreme scale-out workshop: BQCD (Quantum Physics), SeisSol (Geophysics, Seismics), GPI-2/GASPI (Toolkit for HPC), Seven-League Hydro (Astrophysics), ILBDC (Lattice Boltzmann CFD), Iphigenie (Molecular Dynamic), FLASH (Astrophysics), GADGET (Cosmological Dynamics), PSC (Plasma Physics), waLBerla (Lattice Boltzmann CFD), Musubi (Lattice Boltzmann CFD), Vertex3D (Stellar Astrophysics), CIAO (Combustion CFD), and LS1-Mardyn (Material Science). The projects were allowed to use the machine exclusively during the 28 day period, which corresponds to a total of 63.4 million core-hours, of which 43.8 million core-hours were used by the applications, resulting in a utilization of 69%. The top 3 users were using 15.2, 6.4, and 4.7 million core-hours, respectively.

We present, for the first time in the literature, a full reconstruction of the total (linear and non-linear) ISW/Rees-Sciama effect in the presence of massive neutrinos, together with its cross-correlations with CMB-lensing and weak-lensing signals. The present analyses make use of all-sky maps extracted via ray-tracing across the gravitational potential distribution provided by the "Dark Energy and Massive Neutrino Universe" (DEMNUni) project, a set of large-volume, high-resolution cosmological N-body simulations, where neutrinos are treated as separate collisionless particles. We correctly recover, at $1-2\%$ accuracy, the linear predictions from CAMB. Concerning the CMB-lensing and weak-lensing signals, we also recover, with similar accuracy, the signal predicted by Boltzmann codes, once non-linear neutrino corrections to Halofit are accounted for. Interestingly, in the ISW/Rees-Sciama signal, and its cross correlation with lensing, we find an excess of power with respect to the massless case, due to free streaming neutrinos, roughly at the transition scale between the linear and non-linear regimes. The excess is $\sim 5-10\%$ at $l\sim 100$ for the ISW/Rees-Sciama auto power spectrum, depending on the total neutrino mass $M_\nu$, and becomes a factor of $\sim 4$ for $M_\nu=0.3$ eV, at $l\sim 600$, for the ISW/Rees-Sciama cross power with CMB-lensing. This effect should be taken into account for the correct estimation of the CMB temperature bispectrum in the presence of massive neutrinos.

We present radiation transfer models of rotating young stellar objects (YSOs) with hotspots in their atmospheres, inner disk warps and other 3-D effects in the nearby circumstellar environment. Our models are based on the geometry expected from the magneto-accretion theory, where material moving inward in the disk flows along magnetic field lines to the star and creates stellar hotspots upon impact. Due to rotation of the star and magnetosphere, the disk is variably illuminated. We compare our model light curves to data from the Spitzer YSOVAR project (Morales-Calderon et al. 2014, Cody et al. 2014) to determine if these processes can explain the variability observed at optical and mid-infrared wavelengths in young stars. We focus on those variables exhibiting "dipper" behavior that may be periodic, quasi-periodic, or aperiodic. We find that the stellar hotspot size and temperature affects the optical and near-infrared light curves, while the shape and vertical extent of the inner disk warp affects the mid-IR light curve variations. Clumpy disk distributions with non-uniform fractal density structure produce more stochastic light curves. We conclude that the magneto-accretion theory is consistent with certain aspects of the multi-wavelength photometric variability exhibited by low-mass YSOs. More detailed modeling of individual sources can be used to better determine the stellar hotspot and inner disk geometries of particular sources.

We present results from multifrequency radiative hydrodynamical chemistry simulations addressing primordial star formation and related stellar feedback from various populations of stars, stellar energy distributions (SEDs) and initial mass functions. Spectra for massive stars, intermediate-mass stars and regular solar-like stars are adopted over a grid of 150 frequency bins and consistently coupled with hydrodynamics, heavy-element pollution and non-equilibrium species calculations. Powerful massive population III stars are found to be able to largely ionize H and, subsequently, He and He$^+$, causing an inversion of the equation of state and a boost of the Jeans masses in the early intergalactic medium. Radiative effects on star formation rates are between a factor of a few and 1 dex, depending on the SED. Radiative processes are responsible for gas heating and photoevaporation, although emission from soft SEDs has minor impacts. These findings have implications for cosmic gas preheating, primordial direct-collapse black holes, the build-up of "cosmic fossils" such as low-mass dwarf galaxies, the role of AGNi during reionization, the early formation of extended disks and angular-momentum catastrophe.

In this paper we present a large database of weak lensing light cones constructed using different snapshots from the Big MultiDark simulation (BigMDPL). The ray-tracing through different multiple planes has been performed with the GLAMER code accounting both for single source redshifts and for sources distributed along the cosmic time. This first paper presents weak lensing forecasts and results according to the geometry of the VIPERS-W1 and VIPERS-W4 field of view. Additional fields will be available on our database and new ones can be run upon request. Our database also contains some tools for lensing analysis. In this paper we present results for convergence power spectra, one point and high order weak lensing statistics useful for forecasts and for cosmological studies. Covariance matrices have also been computed for the different realisations of the W1 and W4 fields. In addition we compute also galaxy-shear and projected density contrasts for different halo masses at two lens redshifts according to the CFHTLS source redshift distribution both using stacking and cross-correlation techniques, finding very good agreement.

We perform multi-plane ray-tracing using the GLAMER gravitational lensing code within high-resolution light-cones extracted from the CoDECS simulations: a suite of cosmological runs featuring a coupling between Dark Energy and Cold Dark Matter. We show that the presence of the coupling is evident not only in the redshift evolution of the normalisation of the convergence power spectrum, but also in differences in non-linear structure formation with respect to {\Lambda}CDM. Using a tomographic approach under the assumption of a {\Lambda}CDM cosmology, we demonstrate that weak lensing measurements would result in a {\sigma}8 value that changes with the source redshift if the true underlying cosmology is a coupled Dark Energy one. This provides a generic null test for these types of models. We also find that different models of coupled Dark Energy can show either an enhanced or a suppressed correlation between convergence maps with differing source redshifts as compared to {\Lambda}CDM. This would provide a direct way to discriminate between different possible realisations of the coupled Dark Energy scenario. Finally, we discuss the impact of the coupling on several lensing observables for different source redshifts and angular scales with realistic source redshift distributions for current ground-based and future space-based lensing surveys.

We present an algorithm using Principal Component Analysis (PCA) to subtract galaxies from imaging data, and also two algorithms to find strong, galaxy-scale gravitational lenses in the resulting residual image. The combined method is optimized to find full or partial Einstein rings. Starting from a pre-selection of potential massive galaxies, we first perform a PCA to build a set of basis vectors. The galaxy images are reconstructed using the PCA basis and subtracted from the data. We then filter the residual image with two different methods. The first uses a curvelet (curved wavelets) filter of the residual images to enhance any curved/ring feature. The resulting image is transformed in polar coordinates, centered on the lens galaxy center. In these coordinates, a ring is turned into a line, allowing us to detect very faint rings by taking advantage of the integrated signal-to-noise in the ring (a line in polar coordinates). The second way of analysing the PCA-subtracted images identifies structures in the residual images and assesses whether they are lensed images according to their orientation, multiplicity and elongation. We apply the two methods to a sample of simulated Einstein rings, as they would be observed with the ESA Euclid satellite in the VIS band. The polar coordinates transform allows us to reach a completeness of 90% and a purity of 86%, as soon as the signal-to-noise integrated in the ring is higher than 30, and almost independent of the size of the Einstein ring. Finally, we show with real data that our PCA-based galaxy subtraction scheme performs better than traditional subtraction based on model fitting to the data. Our algorithm can be developed and improved further using machine learning and dictionary learning methods, which would extend the capabilities of the method to more complex and diverse galaxy shapes.

The relative impact of radiation pressure and photoionization feedback from young stars on surrounding gas is studied with hydrodynamic radiative transfer (RT) simulations. The calculations focus on the single-scattering (direct radiation pressure) and optically thick regime, and adopt a moment-based RT-method implemented in the moving-mesh code AREPO. The source luminosity, gas density profile and initial temperature are varied. At typical temperatures and densities of molecular clouds, radiation pressure drives velocities of order ~20 km/s over 1-5 Myr; enough to unbind the smaller clouds. However, these estimates ignore the effects of photoionization that naturally occur concurrently. When radiation pressure and photoionization act together, the latter is substantially more efficient, inducing velocities comparable to the sound speed of the hot ionized medium (10-15 km/s) on timescales far shorter than required for accumulating similar momentum with radiation pressure. This mismatch allows photoionization to dominate the feedback as the heating and expansion of gas lowers the central densities, further diminishing the impact of radiation pressure. Our results indicate that a proper treatment of the impact of young stars on the interstellar medium needs to primarily account for their ionization power whereas direct radiation pressure appears to be a secondary effect. This conclusion may change if extreme boosts of the radiation pressure by photon trapping are assumed.

A computer code is described for the simulation of gravitational lensing data. The code incorporates adaptive mesh refinement in choosing which rays to shoot based on the requirements of the source size, location and surface brightness distribution or to find critical curves/caustics. A variety of source surface brightness models are implemented to represent galaxies and quasar emission regions. The lensing mass can be represented by point masses (stars), smoothed simulation particles, analytic halo models, pixelized mass maps or any combination of these. The deflection and beam distortions (convergence and shear) are calculated by modified tree algorithm when halos, point masses or particles are used and by FFT when mass maps are used. The combination of these methods allow for a very large dynamical range to be represented in a single simulation. Individual images of galaxies can be represented in a simulation that covers many square degrees. For an individual strongly lensed quasar, source sizes from the size of the quasar's host galaxy (~ 100 kpc) down to microlensing scales (~ 10^-4 pc) can be probed in a self consistent simulation. Descriptions of various tests of the code's accuracy are given.

We present an extension to multiple planes of the gravitational lensing code {\small GLAMER}. The method entails projecting the mass in the observed light-cone onto a discrete number of lens planes and inverse ray-shooting from the image to the source plane. The mass on each plane can be represented as halos, simulation particles, a projected mass map extracted form a numerical simulation or any combination of these. The image finding is done in a source oriented fashion, where only regions of interest are iteratively refined on an initially coarse image plane grid. The calculations are performed in parallel on shared memory machines. The code is able to handle different types of analytic halos (NFW, NSIE, power-law, etc.), haloes extracted from numerical simulations and clusters constructed from semi-analytic models ({\small MOKA}). Likewise, there are several different options for modeling the source(s) which can be distributed throughout the light-cone. The distribution of matter in the light-cone can be either taken from a pre-existing N-body numerical simulations, from halo catalogs, or are generated from an analytic mass function. We present several tests of the code and demonstrate some of its applications such as generating mock images of galaxy and galaxy cluster lenses.

Mass and spin are the only two parameters needed to completely characterize black holes in General Relativity. However, the interaction between black holes and their environment is where complexity lies, as the relevant physical processes occur over a large range of scales. That is particularly relevant in the case of super-massive black holes (SMBHs), hosted in galaxy centers, and surrounded by swirling gas and various generations of stars. These compete with the SMBH for gas consumption and affect both dynamics and thermodynamics of the gas itself. How the behavior of such fiery environment influence the angular momentum of the gas accreted onto SMBHs, and, hence, black-hole spins is uncertain. We explore the interaction between SMBHs and their environment via first 3D sub-parsec resolution simulations (ranging from 0.1 pc to 1 kpc scales) that study the evolution of the SMBH spin by including the effects of star formation, stellar feedback, radiative transfer, and metal pollution according to the proper stellar yields and lifetimes. This approach is crucial to investigate the impact of star formation processes and feedback effects on the angular momentum of the material that could accrete on the central hole. We find that star formation and feedback mechanisms can locally inject significant amounts of entropy in the surrounding medium, and impact on the inflow inclination angles and Eddington fractions. As a consequence, the resulting trends show upper-intermediate equilibrium values for the spin parameter, a, of about 0.6 - 0.9, corresponding to radiative efficiencies \epsilon = 9% - 15%. These results suggest that star formation feedback taking place in the circum-nuclear disk during the in-fall cannot induce alone very strong chaotic trends in the gas flow, quite independently from the different numerical parameters.

We present results from self-consistent 3D numerical simulations of cosmic structure formation with a multi-frequency radiative transfer scheme and non-equilibrium molecular chemistry of 13 primordial species (e-, H, H+, H-, He, He+, He++, H2, H2+, D, D+, HD, HeH+), performed by using the simulation code GADGET. We describe our implementation and show tests for ionized sphere expansion in a static and dynamic density field around a central radiative source, and for cosmological abundance evolution coupled with the cosmic microwave background radiation. As a demonstrative application of radiative feedback on molecular gas, we run also cosmological simulations of early structure formation in a ~1Mpc size box. Our tests agree well with analytical and numerical expectations. Consistently with other works, we find that ionization fronts from central sources can boost H2 fractions in shock-compressed gas. The tight dependence on H2 lead to a corresponding boost of HD fractions, as well. We see a strong lowering of the the typical molecular abundances up to several orders of magnitudes which partially hinders further gas collapse of pristine neutral gas, and clearly suggests the need of re-ionized gas or metal cooling for the formation of the following generation of structures.

We carry out hydrodynamical simulations of galaxy formation that simultaneously follow radiative transfer of hydrogen-ionizing photons, based on the optically-thin variable Eddington tensor approximation as implemented in the {\small GADGET} code. We consider only star-forming galaxies as sources and examine to what extent they can yield a reasonable reionization history and thermal state of the intergalactic medium at redshifts around $z\sim 3$. This serves as an important benchmark for our self-consistent methodology to simulate galaxy formation and reionization, and for future improvements through accounting of other sources and other wavelength ranges. We find that star formation alone is sufficient for reionizing the Universe by redshift $z\sim6$. For a suitable choice of the escape fraction and the heating efficiency, our models are approximately able to account at the same time for the one-point function and the power spectrum of the Lyman-$\alpha$ forest. The radiation field has an important impact on the star formation rate density in our simulations and significantly lowers the gaseous and stellar fractions in low-mass dark matter halos. Our results thus directly demonstrate the importance of radiative feedback for galaxy formation. The spatial and temporal importance of this effect can be studied accurately with the modelling technique explored here, allowing more faithful simulations of galaxy formation.

We present a numerical implementation of radiative transfer based on an explicitly photon-conserving advection scheme, where radiative fluxes over the cell interfaces of a structured or unstructured mesh are calculated with a second-order reconstruction of the intensity field. The approach employs a direct discretisation of the radiative transfer equation in Boltzmann form with adjustable angular resolution that in principle works equally well in the optically thin and optically thick regimes. In our most general formulation of the scheme, the local radiation field is decomposed into a linear sum of directional bins of equal solid-angle, tessellating the unit sphere. Each of these "cone-fields" is transported independently, with constant intensity as a function of direction within the cone. Photons propagate at the speed of light (or optionally using a reduced speed of light approximation to allow larger timesteps), yielding a fully time-dependent solution of the radiative transfer equation that can naturally cope with an arbitrary number of sources, as well as with scattering. The method casts sharp shadows, subject to the limitations induced by the adopted angular resolution. If the number of point sources is small and scattering is unimportant, our implementation can alternatively treat each source exactly in angular space, producing shadows whose sharpness is only limited by the grid resolution. A third hybrid alternative is to treat only a small number of the locally most luminous point sources explicitly, with the rest of the radiation intensity followed in a radiative diffusion approximation. We have implemented the method in the moving-mesh code {\small AREPO}, where it is coupled to the hydrodynamics in an operator splitting approach that subcycles the radiative transfer alternatingly with the hydrodynamical evolution steps.

We present a novel numerical implementation of radiative transfer in the cosmological smoothed particle hydrodynamics (SPH) simulation code {\small GADGET}. It is based on a fast, robust and photon-conserving integration scheme where the radiation transport problem is approximated in terms of moments of the transfer equation and by using a variable Eddington tensor as a closure relation, following the `OTVET'-suggestion of Gnedin & Abel. We derive a suitable anisotropic diffusion operator for use in the SPH discretization of the local photon transport, and we combine this with an implicit solver that guarantees robustness and photon conservation. This entails a matrix inversion problem of a huge, sparsely populated matrix that is distributed in memory in our parallel code. We solve this task iteratively with a conjugate gradient scheme. Finally, to model photon sink processes we consider ionisation and recombination processes of hydrogen, which is represented with a chemical network that is evolved with an implicit time integration scheme. We present several tests of our implementation, including single and multiple sources in static uniform density fields with and without temperature evolution, shadowing by a dense clump, and multiple sources in a static cosmological density field. All tests agree quite well with analytical computations or with predictions from other radiative transfer codes, except for shadowing. However, unlike most other radiative transfer codes presently in use for studying reionisation, our new method can be used on-the-fly during dynamical cosmological simulation, allowing simultaneous treatments of galaxy formation and the reionisation process of the Universe.