© The Authors 2023

1. Introduction

Feedback models for galaxy evolution are mostly based on theoretical predictions, hydrodynamical simulations, and kinematical models based on observations (Fabian 2012; Morganti 2017; Kudritzki et al. 2021). Therefore, an in-depth comprehension of multi-phase gas kinematics in galaxies is crucial in order to determine the mutual interaction between the phases of interstellar medium (ISM) and to shed light on energy injection and ejection mechanisms. Integral Field Unit (IFU) observations, together with kinematical models, allow us to constrain with unprecedented detail and accuracy the 2D and 3D kinematics of gas discs in galaxies (Oh et al. 2008; Di Teodoro & Fraternali 2015; Bouché et al. 2015). Despite considerable advances in the modelling of gas kinematics in galaxy discs, galaxy-wide multi-phase outflows are still poorly constrained because of their complex observed geometry, kinematics and surface brightness distribution (e.g., Veilleux et al. 2020).

The unified model of active galactic nuclei (AGN) predicts that radiation from the accretion disc around the supermassive black hole is collimated by a torus of obscuring dust and gas, which is broadly symmetric around the accretion flow axis (Urry & Padovani 1995). This radiation impinges on the ISM, ionising the gas and transferring mechanical energy through radiation pressure and disc winds. Consequently, the outflowing gas extends outward from the nucleus assuming a typical (bi)conical geometry. These cones are more frequently observed in local Seyfert galaxies, and more often in type-II (obscured) than type-I (unobscured) Seyferts, which is consistent with the unification model (Schmitt & Kinney 2016). The best tracer of the ionised phase of these outflows on 10² − 10³ pc scales is the forbidden emission line doublet [OIII]λλ4959, 5007, as it can only be produced in low-density regions and therefore does not trace the subparsec scales of the broad line region (BLR). In the presence of outflows, the spectral profile of the [OIII]λ5007 emission line typically assumes an asymmetric shape, with a broad, blueshifted wing¹ (e.g., Boroson 2005; Zakamska & Greene 2014; Perna et al. 2017; Bae & Woo 2018). A detailed study of gas kinematics and ionisation mechanisms in these cones is extremely important in order to shed light on the interaction between AGN and their host galaxies, and on the multi-phase nature of outflows (e.g., Crenshaw et al. 2015; Karouzos et al. 2016; Bae & Woo 2018; Cresci & Maiolino 2018; Venturi et al. 2018; Mingozzi et al. 2019; Fluetsch et al. 2019; Kraemer et al. 2020; Marasco et al. 2020; Tozzi et al. 2021). AGN-driven outflows represent the main manifestation of radiative feedback in action, as they have the potential to sweep away the gas content of the host galaxy, thus quenching star formation inside the outflow cavity from low to high redshift (e.g., Feruglio et al. 2010; Hopkins & Elvis 2010; Sturm et al. 2011; Cano-Díaz et al. 2012; Zubovas & King 2012; Cicone et al. 2014; Carniani et al. 2015, 2016; Cresci et al. 2015a). As a result, they can shape galaxies and cause the well-known empirical relations between black hole mass (M_BH) and host galaxy properties, that is, the mass or luminosity of the bulge (M_bulge and L_bulge, respectively) and velocity dispersion (σ) (Gebhardt et al. 2000; Marconi & Hunt 2003; Ferrarese & Ford 2005).

Even though AGN-powered outflows are thought to play a key role in galaxy evolution, an accurate determination of outflow properties and driving mechanisms is still missing, which is due to the large uncertainties related with the estimated outflow properties in different gas phases (Cicone et al. 2014; Veilleux et al. 2017; Harrison et al. 2018). In most sources, because of unresolved observations, it is unclear whether the outflow morphology is conical or shell-like, and the physical processes that drive the coupling of the energy and momentum released by the AGN with the surrounding ISM are a matter of significant debate. The most accepted model predicts that a highly ionised nuclear wind accelerated by the BH radiation pressure shocks the ISM of the host galaxy, creating an expanding bubble of hot gas, with the post-shock medium going through a momentum-conserving phase followed by an adiabatically expanding phase (King 2003; King et al. 2011; Cicone et al. 2014; Feruglio et al. 2015; Tombesi et al. 2015). To address these debated topics, many kinematical models have been proposed to analyse observations of galactic outflows, but so far only a few have been tested with spatially resolved integral field spectroscopy (IFS) data (e.g., the AGN-outflow model of VLT/SINFONI data from Müller-Sánchez et al. 2011). Simpler models based on long-slit observations have been proposed to investigate the kinematics of the narrow-line region (NLR) in local Seyfert galaxies, modelling the position-velocity (PV) diagram (Crenshaw et al. 2000; Crenshaw & Kraemer 2000; Das et al. 2005; Fischer et al. 2010; Storchi-Bergmann et al. 2010).

Venturi et al. (2017) made a first comparison between the observed moment maps obtained from MUSE data of NGC 4945 and NGC 1365, and a kinematic toy model to investigate the outflow kinematics. This toy model is the starting point for our current kinematic model. Up to now, literature models, such as those mentioned above, have been computed by adopting a conical outflow geometry, with velocity field and gas surface brightness of the emitting gas parameterised by analytical functions of the distance from the cone vertex. The velocity profile, at a given position on the sky is given by:

(1)

where P is a given position on the sky, and Σ_P(s, v) is the corresponding gas surface brightness distribution at a given velocity v and coordinate s along the line of sight (LOS). The best parameters are then found to reproduce f_P(v) for any given spaxel, or its average velocity and velocity dispersion. In a few cases, the model takes into account the finite spatial and spectral resolution by convolving the data with the proper smoothing kernels (e.g., Bae & Woo 2016; Karouzos et al. 2016). However, the assumption of a smooth function Σ(s, v) to describe the gas emissivity is in contrast with observations, which in nearby galaxies show clumpy and irregular structures. Therefore, an important question to be addressed is whether the complex velocity field observed in nearby sources, which has never been reproduced by previous studies, is the result of a complex velocity distribution within the propagating outflow, or the consequence of a clumpy distribution of ionised gas clouds.

To address these debated questions, we propose an innovative approach to modelling gas kinematics and determining outflow properties, explaining improvements with respect to previous methods. In this paper, we present our new kinematic model MOKA^3D (Modelling Outflows and Kinematics of Agn in 3D), discuss its main degeneracies and show examples of its application to three Seyfert II galaxies from the MAGNUM survey (Measuring Active Galactic Nuclei Under MUSE Microscope, Cresci et al. 2015b; Venturi et al. 2018, 2021; Mingozzi et al. 2019), NGC 4945, Circinus, and NGC 7582.

The paper is constructed as follows. In Sect. 2, we describe step by step the 3D model operation. In Sect. 3 we present the application of our method to simulated data. In Sect. 4, we introduce the sample to which we applied the model, describing the gas kinematic features, and the results for the kinematics and energetics obtained following the model application. In Sect. 5, we summarise our results and present future developments.

2. Kinematic model

In this section, we present an overview of our model, describing the main steps to infer the kinematics and orientation of observed galactic structures. In Fig. 1, we show a schematic flowchart of the MOKA^3D model operation and fitting, which are described in the following sections.

Fig. 1. Flowchart of our MOKA3D kinematic model. The method consists in assigning the model parameters and assuming a coordinate system and a velocity field. Each cloud is assigned the coordinates in the source reference frame (green), which is later transformed to the observer reference frame (blue), through Euler rotation matrices. If necessary, the method requires the convolution with the PSF and/or the spectral LSF. Then, the model is binned in the space sampled by the observed data cube (red). The model emission line profile percentiles at 1% and 99% are compared to the observed ones and, if they are comparable, each particle is weighted according to the observed emission in the corresponding bin. The model is evaluated by means of the κ estimator (Eq. (6)). The free parameters are varied and the loop repeats, until finding the suitable parameter configuration that minimises κ and reproduces the observed spectrum in each spaxel.

2.1. Model setup

We adopt a spherical coordinate system and create a uniform distribution of 10⁷ point-like synthetic emitting sources (clouds, hereafter) distributed according to an input geometrical distribution chosen by the user (e.g., a cone when modeling outflows). We decided to adopt this number of clouds after considering the average MUSE emission line data cube size, that is ∼300 × 300 × 40. In particular, as discussed in Sect. 2.3, we need at least a few model clouds per voxel (volume pixel) to properly model the observed features. The spherical coordinate system is selected to model radial outflows, but the method works with any kind of reference system, such as cylindrical or Cartesian. The 3D spatial location of each cloud is specified by three coordinates: the semi-polar angle θ measured from the cone axis (0° ≤θ ≤ 180°); the azimuthal angle ϕ of the cloud orthogonal projection on a reference plane passing through the origin and orthogonal to the axis (measured clockwise, 0° ≤ϕ ≤ 360°); and the radial distance r from the origin (in arcsec). We start by assuming the surface brightness distribution of the emitting clouds as follows:

(2)

where f₀ is the flux value at the apex of the cone and r₀ is an arbitrary scale-radius. For what concerns the outflow velocity, we start by assuming it as constant with radius (V(r) = V₀).

2.2. Reference frame

We define the source reference system in Cartesian coordinates as (X, Y, Z) = (sin(θ)cos(ϕ),sin(θ)sin(ϕ),cos(θ)). Therefore, the status of each cloud is specified by position and velocity vectors, namely, P = (X, Y, Z) and V = (V_X, V_Y, V_Z), with the velocity vector defined as:

(3)

It is also possible to associate a random velocity dispersion component σ_rand to each cloud, but in this work we have not taken advantage of this possibility. Then, we consider the observer’s frame (x, y, z), where (x, y) are the coordinates on the plane of the sky and z is directed along the LOS. The source and observer reference frames are shown in Fig. 2, as blue and red, respectively.

Fig. 2. Source (blue) and observer (red) reference frame. Euler angles α, β, and γ are used to transform from the source frame to the observer’s frame and vice versa. The purple dotted line represents the line of nodes, N. The observer z axis coincides with the LOS, and (x, y) is the plane of the sky.

We transform P and V from the source to the reference frame by means of Euler rotation matrices:

(4)

where R_γ, R_β, R_α are the Euler rotation matrices, with the rotation angles shown in Fig. 2 and defined as follows:

• α: corresponds to the angular separation between the line of nodes and the right ascension coordinate in the observer reference system. The α angle is irrelevant if the sobserved source is axially symmetric around the source z-axis.

• β: corresponds to the outflow axis inclination with respect to the plane of the sky (e.g., β = 180° in the case of a source pointing away from the observer, with the z-axis along the LOS; β = 90° for a source axis lying on the plane of the sky).

• γ: corresponds to the projected source major axis inclination in the plane of the sky with respect to the line of nodes, indicating the rotation direction and is also known as the position angle (P.A.).

The velocity component along the LOS in the observer reference frame is defined as:

(5)

where n is the unit vector pointing away from the observer, and v_sys is the outflow systemic velocity with respect to the galaxy. Once the clouds are projected, the model allows both a spatial and spectral convolution to account for observational effects. The observed position (x, y) of each cloud is randomly shifted on the plane of the sky according to a 2D probability distribution given by the PSF measured from the data. If necessary, to account for the intrinsic spectral resolution, the observed LOS velocity can be similarly convolved with the line spread function (LSF), but in this work we have not taken advantage of this possibility. In the following, we refer to the outcome of this step as the unweighted model.

2.3. Unweighted model

At this stage, the unweighted model is binned in the observed (x, y, v_z) space, which is the one sampled in data cubes produced by IFU observations. The process is schematised in Fig. 3. The spatial and spectral extensions of the model cube are forced to match the spatial and spectral sampling of the data, respectively. A voxel of the obtained 3D model cube, with two spatial extensions and a spectral-velocity one, uniformly populated with clouds, is highlighted in black in Fig. 3. In this phase, all clouds in the same bin have a weight associated to our initially assumed radial-dependent flux function (Eq. (2)). In particular, all the clouds corresponding to the same spaxel have the same intensity, being dependent only on the radial distance from the vertex of the cone. Spreading 10⁷ clouds homogeneously in the binned model cube guarantees that at least a few clouds populate each voxel, and therefore that the weighting procedure described in the following Sect. 2.4 succeeds. Since multiple clouds can fall in the same bin, as shown in Fig. 3, their contribution to the emission of the corresponding model spectral channel is assumed to be equal.

Fig. 3. Coordinate transformation and binning procedure of the model: from the simple conical geometry in the source reference frame, uniformly populated with clouds, to the 3D model cube in the observer’s reference frame. The model is then binned to assign to each cloud a position in the (x, y, vz) space, as highlighted for the purple cloud. The model cube has the same spatial and spectral resolution as the observed data cube in order to weight the clouds according to the observed flux in the corresponding spaxel. Dashed grey lines show two adjacent voxels, corresponding to different spatial and/or spectral channels.

Once the unweighted model is computed, each spaxel (spectral pixel) has an associated spectrum (see Eq. (1)), which depends on the model geometry and the parameters of the flux and velocity functions we assumed. This model spectrum can be compared with the corresponding observed spectrum in the same spaxel. It is then possible to compute projected flux, velocity, and velocity dispersion maps from the momenta of the model cube. Figure 4 shows two examples of unweighted models: a filled bi-conical outflow (top panels) and a hollow conical one (bottom panels).

Fig. 4. Example of unweighted models simulating an outflow with the spatial and spectral resolution of a galaxy observed with MUSE at z = 0.001. From left to right: Line-velocity integrated map, LOS velocity map, and velocity dispersion map. Top panels: full biconical model with P.A. = 40°, inclination with respect to the plane of the sky of β = 85°, a constant radial velocity of 1200 km s−1, and a cone semi-aperture of θOUT = 30°. Bottom panels: single hollow conical model with P.A. = 90°, β = 65°, and inner and outer semi-aperture θIN = 25° ,θOUT = 35°. Both models have been convolved with a spatial PSF = 0.3″ and no spectral convolution.

2.4. Weighted model

We use the line flux F_obs from each (x, y, v_z) voxel of the observed data cube, to determine the surface brightness of each model cloud. In particular, we assign a weight w = F_obs/N_mod to each cloud within the volume element specified by the (x, y, v_z) coordinates, where N_mod is the total number of clouds ending up in that voxel. We want to clarify that, hereafter, each cloud is assigned a weight that is independent of the distance from the vertex of the cone. Therefore, each cloud emission is not related to previous assumptions on the analytical surface brightness distribution (this distribution could be any function; see e.g., Eq. (2)). As an example, all the clouds in the highlighted black voxel in Fig. 3 have the same weight because they belong to the same voxel, and are therefore weighted with the same flux. We can now compare the observed line profile in each spaxel with those obtained from a variety of weighted model cubes, each one computed using a different combination of geometrical and physical parameters. Assigning a weight to each cloud ensures that the model velocity profile matches the observed one by construction, provided that model clouds occupy the same velocity range spanned by the data. The latter condition will not be satisfied if the model geometry and velocity field is incorrect, which gives us the means to constrain both the cone geometry and the intrinsic velocity field at the same time. The procedure to infer the outflow parameters through a fitting procedure is described in the following section.

2.5. Fit of the model to observed data

The fitting algorithm consists in a loop over the parameter space until the best set of parameters is obtained, which is that which minimises the difference between the observed and modelled emission in each voxel. First, we define the free parameters we want to explore for the MUSE example: intrinsic radial outflow velocity, outflow systemic velocity with respect to the galaxy systemic velocity (V and v_sys in Eq. (5), respectively), and inclination with respect to the LOS. Then, following the procedure shown in Fig. 1, for each set of parameters we create a weighted model cube, whose emission is compared to the observed emission over the velocity range defined by the 1% and 99% percentiles of the observed LOS velocity distribution, v₁ and v₉₉, respectively. Finally, we evaluate how well each model reproduces the observed spectrum by means of a customised goodness of fit estimator, defined as:

(6)

Here, v_i represents the ith spectral channel of the jth spaxel. Therefore, S_Oj(v_i) and S_Mj(v_i) are the observed and model spectral flux in the ith voxel, respectively, and δs_j is the uncertainty on the flux, which is assumed to be constant in all spectral channels. We minimise κ to find the best set of model parameters. Although the κ estimator is defined as the standard χ² estimator, its definition relies on different assumptions. For this reason, we decided to refer to it as κ instead of χ². A proper computation of the free parameter uncertainties and the development of a tailored statistic will be addressed in future work.

Figure 5 shows a comparison between the observed (black) [NII]λ6584 emission line profile extracted from a spaxel with a prominent blueshifted wing in NGC 4945 and the model (red). The top panels show the unweighted model emission for the best-fit model (left) and a model with incorrect inclination (right). The bottom panels show the weighted version of the top panels. The best model configuration on the bottom left is obtained with the parameters listed in Table 1, and perfectly reproduces the observed emission in each velocity channel. In the right panels, a different set of model parameters is not able to reproduce the observed blue wing, for either the weighted or unweighted model. This is due to the fact that, for this incorrect set of parameters, no cloud contributes to observed LOS velocities lower than −300 km s⁻¹. Indeed, for a model to be successful, it is crucial that all spaxels contain model clouds populating all velocity channels where emission is observed.

Fig. 5. Comparison between the observed [NII]λ6584 emission line profile (black) and the model emission (red) of NGC 4945, extracted from a spaxel with a prominent blue-shifted wing. Top and bottom panels show the unweighted and weighted emission, respectively. Left panels: the model emission is obtained with the best fit parameters reported in Table 1. Right panels: as a demonstrative example, we show the model line profile as obtained by assuming an unsuitable inclination and opening angle. Such a configuration is clearly not able to to reproduce the high-velocity, blue-shifted emission, thus causing a sharp cut off. The absence of clouds in all velocity channels where emission is observed, prevent the model to reproduce the observed outflow properties.

2.6. Solving the degeneracies

The main challenge when creating a 3D kinematical model is to recognise the degeneracies affecting the fit results. The main degeneracy affecting MOKA^3D is that very high outflow velocities, combined with a range of different inclinations towards the observer’s LOS, always allow us to reproduce the observed line profile in each spaxel, even with an incorrect set of geometrical parameters, because all the velocity channels of each spaxel will always be populated with model clouds.

As an example of this, Fig. 6 shows the moment maps of observed data (top panels), and two models with a radial outflow velocity of 3000 km s⁻¹, which is above the observed maximum velocity of ∼900 km s⁻¹. The model in the middle panels is obtained with β = 50°, v_sys = +1000 km s⁻¹, and an outer semi-opening angle of θ_out = 48°; the model in the bottom panels instead is obtained with β = 120°, v_sys = −500 km s⁻¹, and an outer semi-opening angle of θ_out = 38°.

Fig. 6. Comparison of observed data and two degenerate models of NGC 4945. Moment maps of observed data (top panels) and of two models with high radial outflow velocity and incorrect inclination with respect to the plane of the sky, incorrect outer opening angle, and incorrect systemic velocity (middle and bottom panels) with respect to the best model shown in Fig. 12. See Sect. 2.6 for details of the parameters of the degenerate models and how to overcome them.

Figure 6 shows that all models with high radial outflow velocity are able to reproduce the observed emission and kinematics in each spaxel equally well, even with a very different geometrical configuration. If the outflow cone aperture is wide enough to reproduce the observed cone aperture in the plane of the sky, the combination of high velocities and a wide range of inclination angles with respect to the plane of the sky can provide a wide range of velocities along the LOS, fitting any kind of observed data once the weighting scheme described above is applied.

As each cloud is assigned a weight, which is measured from the observed data and is dependent on the (x, y, v_z) bin where the cloud falls, there are two configurations causing a cloud to be assigned a zero weight: either the cloud has a velocity v_z beyond the observed velocity ranges, or no emission is detected in the cloud 3D position-velocity. Therefore, an outflow radial velocity that results in velocity percentiles well above the observed boundaries v₁ and v₉₉, as in the cases shown in Fig. 6, allows the model to reproduce the observed spectrum in each spaxel by assigning zero weight to clouds at the model edges, producing a flattened model on the plane of the sky.

Figure 7 schematically shows the effect on the distribution of model clouds by creating a model with much higher radial outflow velocity than that observed. The solid and dotted black lines define the intrinsic outflow aperture and the LOS direction, respectively. The grey dotted lines mark the region containing clouds with projected velocity within the v₁, v₉₉ boundaries. The grey clouds, as opposed to the green ones, have projected velocities above the observed boundaries, and therefore the model assigns them a null weight. The top and bottom panels show the side and top view of the model, respectively. The larger the outflow velocity with respect to the maximum observed velocities, the more flattened the distribution of model clouds with non-zero weight on the plane of the sky.

Fig. 7. Schematic explanation of the effect of an excessive radial outflow velocity with fixed model aperture and inclination. The solid black lines define the outflowing cone aperture, the dotted black lines show the LOS directed towards the observer, and the grey dotted lines show the region containing the clouds with projected velocity within the v1, v99 observed boundaries. The grey clouds have LOS velocity above v1, v99, i.e. their weights are zero. Black arrows show the decomposition of the outflow velocity. The top panel shows the view from the top of the model, while the bottom panel shows the flattening effect as seen from the plane of the sky.

To overcome this degeneracy issue, we constrained the model parameters as follows. While following the loop over the free parameters (see scheme in Fig. 1), the algorithm discards all the combinations of model parameters that result in having percentile velocities v₁, v₉₉ –derived from the integrated unweighted model spectrum–that are different from the observed ones. This is done by defining a parameter f and imposing for each model that and . Here, f ∈ (0.2,0.05) is a refinement parameter that progressively decreases until reaching f = 0.05. In this way, it is guaranteed that the difference between the observed and model percentile velocities differs by no more than 5%. This constraint then requires a suitable combination of outflow velocity, cone aperture, and inclination, to reproduce the observations, in addition to representing an optimal method to remove the degeneracies.

3. MOKA^3D test on simulated data

In this section, we present the results of the application of MOKA^3D to simulated data cubes in order to test the capabilities and limits of the model. We generated four different simulated cubes with surface brightness distribution and kinematics similar to those observed in nearby and high-redshift sources. The purpose is to simulate the observed spatial and spectral resolution of new generation IFU and to test the degree of reliability of the model with increasing complexity. We created each simulated cube using a spatial and spectral sampling very similar to MUSE data cubes, that is 0.2″/pixel and 55 km s⁻¹, respectively.

We tested different inclinations with respect to the LOS, different inner and outer cone opening angles, and different velocity fields and PSFs in order to simulate the outflow features observed with MUSE in our sample. Adopting different PSF sizes allows us to simulate different instrument resolutions and source distances in order to test the performance of MOKA^3D for both low- and high-z sources. For the surface brightness emission, we adopted a random distribution of 100 clumps, which is a sufficient number to produce moment maps similar to those of the MUSE sample. The clump is a 3D normal distribution of increased flux centred around a randomly selected (r, θ, ϕ) coordinate, between the minimum and maximum of the chosen interval of spherical coordinates. The outflow is characterised by a constant surface brightness distribution, which is made irregular by the presence of the clumps. We want to stress that the number of simulated clumps is irrelevant to the result of the fit.

3.1. Fit of simulated data

For each simulated cube we inferred the best geometrical and kinematic parameters by running the fitting algorithm outlined in Sect. 2, progressively increasing the number of free parameters. The comparison between simulated and best-fit moment maps for the first three models is shown in Fig. 8, and the intrinsic and best-fit parameters are reported in Table 2. Both the best-fit moment maps (Fig. 8) and model parameters (Table 2) are in good agreement with the simulated data. The first moment (velocity) maps of each simulated cube clearly show that assuming a very simple velocity field and irregular flux emission, results in very complex and irregular kinematic features, despite the intrinsic velocity field being regular. For the first three models, the fitted parameters are listed in Table 2; for Model 4, the chosen parameters are reported in Table 3. The parameters labelled with an asterisk were kept fixed during the fit.

Fig. 8. Moment maps and integrated emission comparison between three simulated and best-fit models. From top to bottom: Model 1, Model 2, and Model 3, whose parameters are listed in Table 2. Left panel: comparison of simulated and best-fit moment maps. From left to right: integrated emission, LOS velocity, and velocity dispersion maps. Right panel: comparison of simulated data (dotted blue) and best-fit integrated profile (red). Residual emission is shown in green.

In the right panel of Fig. 8 we show the comparison between simulated and best-fit integrated emission line profiles, plotted in dashed blue and solid red, respectively. The line profile residuals shown with the solid green line are to be mainly ascribed to the discrepancies in model inclination and inner/outer opening angle, which result in missing clouds at the correct 3D position, as explained in Sect. 2.5.

Model 4 was designed to be as more similar as possible to the configuration that we expect in our MAGNUM sample. In particular, we simulated a constant radial velocity outflow characterised by clumpy emission superimposed on a rotating gas disc extending above the field of view (FOV). We tested the MOKA^3D capabilities by fitting the total simulated cube emission with a simple conical model and adopting a constant radial velocity field in order to test whether or not the model is able to derive the real kinematical and geometrical outflow parameters despite the presence of an underlying disc. We simulated a uniform disc emission with an intensity ∼10⁻³ times smaller than the emission of the outflow clumps, and adopted a constant rotating disc velocity of 150 km s⁻¹, which is consistent with the observed stellar and gas rotating velocities in the inner kiloparsec of galaxy discs (Sofue 2015; Yoon et al. 2021; Vukcevic 2021). The disc axis is tilted by 5° with respect to the outflow axis. For Model 4, the total emission was convolved with a spatial PSF of 0.6″.

Figure 9 shows the comparison between simulated and best-fit model moment maps and emission line residuals. The best-fit model in Fig. 9 is obtained by fitting the outflow inclination (β), radial velocity (v_r), inner and outer opening angle (θ_IN and θ_OUT; see Table 3). From this test case it emerges that, in the case of combined emission of multiple kinematical components (disc and outflow), fitting only the dominant outflow component still provides accurate results. The line profile residuals in Fig. 9 show that the best-fit model is missing clouds centred at 0 km s⁻¹, which is due to an underlying rotating disc that was not included in the model.

Fig. 9. Same as Fig. 8 but for Model 4. Left panels: moment maps comparison of simulated data (top) and best-fit model (bottom) obtained with four free parameters. From left to right: integrated emission, LOS velocity, and velocity dispersion maps. Right panel: integrated emission line comparison between simulated (dotted blue) and best-fit model (red), with residuals in green.

We tested models with an extremely wide range of axis inclinations, inner and outer opening angles, random clump distributions, and velocity fields in order to account for a variety of possible configurations. From all these tests, it emerges that MOKA^3D correctly derives the kinematic and geometric parameters we used to create the simulated data cubes, with the typical PSF in low- and high-z data. In particular, our model correctly derives, with unprecedented accuracy, the outflow inclination with respect to the LOS, the 3D velocity field, and the outer opening angle. These represent the fundamental parameters to be measured to achieve accurate estimates of the outflow energetics. As expected, increasing the number of free parameters of the fit causes an increase in the degeneracies, with a corresponding lower accuracy even with the method outlined in Sect. 2.6.

Fitting all the outflow parameters listed in Tables 2 or 3, including those labelled with the asterisk, allows for many possible configurations and unreliable best-fit parameters. For example, for Model 1 and 2, we observe that the major consequence of keeping all the parameters is that a hollow conical best-fit model is obtained, despite the intrinsic inner opening angle being exactly zero. We observed that for full conical simulated data, fixing the inner opening angle to zero results in an optimal estimate of the other parameters (see parameters of Model 1). Setting the inner opening angle as a free parameter instead leads to an incorrect hollow cone configuration (see Model 2). As shown for Model 3 instead, starting from hollow conical simulated data results in a correct estimate of the inner opening angle, without the need to fix any parameter.

Model 4 represents the most reliable test because of the co-presence of disc and outflow emission. For this case, we decided to model only the outflow emission, exactly as we did for the MUSE sample (see Sect. 4). As shown by the best-fit parameters (Table 3) and the moment maps comparison (Fig. 9), MOKA^3D is a reliable tool for deriving the outflow properties with great accuracy.

3.2. MOKA^3D limits

Based on the tests we run on simulated data, it emerges that our MOKA^3D method has limitations in the derivation of the kinematical and geometrical properties of simulated complex outflows data. In particular, based on the results obtained for Model 4 (see Fig. 9 and Table 3), which is created in order to simulate a typical data cube from the MUSE sample, it emerges that degeneracies increase with an increased number of free model parameters. For Model 4, we tested the accuracy of our method by progressively increasing the number of free parameters and computed, voxel by voxel, the difference between the best-fit model and simulated cube. We define the accuracy as:

(7)

where I_{simulated, i} and I_fit, i are the intensity of the ith voxel of the masked simulated cube and the best-fit cube, respectively. N is the total number of unmasked voxels. By definition, with the accuracy A approaching unity, the corresponding model can be considered a good representation of the simulated data cube. Moreover, as a consequence of the weighting procedure (see Sect. 2.4), the intensity of any voxel of the weighted model cube will always be smaller than or equal to the corresponding intensity of the simulated/observed data cube (I_{simulated, i} ≥ I_fit, i). We fitted the simulated cube with an increasing number of free parameters, starting with the radial outflow velocity (v_r), and then adding the inclination with respect to the plane of the sky (β), inner opening angle (θ_IN), outer opening angle (θ_OUT), position angle (P.A.), rotation, and expansion or contraction velocity (v_θ, v_ϕ). For the test case of Model 4, we observed that the accuracy level remains > 95% fitting up to the first four outflow parameters. Adding the fifth free parameter, the accuracy drops to ∼80%. The worst scenario is obtained when fitting all seven parameters, which results in an accuracy level of ∼60%.

The test case of Model 4 is highly representative of what happens when only the outflow conical model is fitted, despite the presence of an underlying rotating galaxy disc. In particular, fitting up to four outflow kinematical and geometrical parameters results in a best-fit model cube that, in each voxel, is in very good agreement with the simulated cube.

Moreover, we observed that MOKA^3D correctly derives up to four outflow parameters – despite the presence of an underlying disc – when the observed outflow to disc flux ratio is ≥10^2.5. Figure 10 shows a comparison of the integrated [OIII]λ5007 emission extracted from a clump of the outflow and disc in the Circinus galaxy. The disc emission is multiplied by 10^2.5. Therefore, adopting in our simulations an average flux ratio of 10^2.5–10³ agrees with the value observed in our MUSE sample. We expect that in cases with lower values of outflow to disc flux ratio, MOKA^3D would not be able to properly recover the intrinsic outflow properties. The combination of outflow and disc fitting will be addressed in future work.

Fig. 10. Example of outflow- and disc-selected clumps in the Circinus galaxy. Top panel: 2D [OIII]λ5007 emission map of the Circinus galaxy. Orange and black circles show the average regions where the disc- and outflow-integrated emission are extracted, respectively. Bottom panel: comparison of outflow and disc-integrated spectra (×102.5) extracted from the corresponding regions shown in the left panel.

4. MOKA^3D application to MAGNUM galaxies

To show the application of MOKA^3D to real sources, we selected a sample of three nearby Seyfert-II galaxies from the MAGNUM survey, namely NGC 4945, Circinus and NGC 7582. In this section, we briefly introduce the main properties and gas kinematical features of this sample.

4.1. Preliminary spectroscopic analysis

Data cubes were analysed by means of a set of custom Python scripts to first subtract stellar continuum, and then fit multiple Gaussian components to the emission lines, thus finally obtaining an emission-line model cube for each emission line. For a more detailed description of the data reduction and the spectroscopic analysis we refer to Mingozzi et al. (2019), Marasco et al. (2020), and Tozzi et al. (2021).

4.2. Ionised emission

For each galaxy, we extracted a subcube centred on the nuclear region where the ionisation cone is observed. In the case of Circinus and NGC 7582, we used the [OIII]λ5007 emission line to map outflow properties; for NGC 4945 instead, we used [NII]λ6584, because [OIII] emission is highly obscured by dust in this edge-on galaxy (see e.g., Fig. A1 in Mingozzi et al. 2019). We computed the S/N of the maps by considering the ratio between the peak of the fitted emission line and the rms of the fit residuals; we then applied a S/N threshold of 3. The subcube was extracted from the emission-line model cube –which corresponds to the sum of all Gaussian components used to fit the line– after deconvolving the [OIII] line profile by the MUSE instrumental broadening². Once we had limited the area we are interested in, we created moment maps of the ionised emission from the cube fit: the integrated line flux (moment of order 0), the flux-weighted LOS velocity (moment of order 1) and the velocity dispersion maps (moment of order 2). Below, we show how this maps are used to provide a first comparison with the final model results. In this way, we want to demonstrate that correctly reproducing the three moment maps is necessary but not sufficient to guarantee that a model is a faithful representation of the outflow features and properties; a more detailed discussion about degeneracies is carried out in Sect. 2.6.

4.3. Moment maps

In Fig. 11, we show the moment maps for each source: starting from the left, integrated emission line, LOS velocity, and velocity dispersion maps; from top to bottom: NGC 4945, Circinus, and NGC 7582.

Fig. 11. Moment maps of the sample. From top to bottom: NGC 4945, Circinus, and NGC 7582. From left to right: Integrated surface brightness map in units of 10−20 erg cm−2 s−1 arcsec−2, intensity-weighted velocity and velocity dispersion. The first row shows the [NII]λ6584 emission in NGC 4945 and the last two rows refers to the [OIII]λ5007 emission. North is up. The maps are 1 spaxel−σ spatially re-smoothed, so as to get a better visual output, and pixels with S/N ≥ 3 are masked. The maps are obtained from the spectroscopic analysis of the modelled emission line from MUSE data (e.g., Marasco et al. 2020).

The selected galaxies were chosen to show extremely clumpy, spatially resolved ionised outflows, as highlighted by the integrated emission maps in Fig. 11. The projected distance covered by the outflow ranges from ∼1 kpc in Circinus to ∼3 kpc in NGC 7582 (Juneau et al. 2022). The intensity maps show a bi-conical axis-symmetric geometry in all cases except for Circinus, where dust and gas in the galaxy disc are probably obscuring the counter cone. The receding cone in NGC 7582 and NGC 4945 is more dust obscured than the approaching one, as is usually observed in Seyfert galaxies and AGN in general. A common feature that we detected in all sources is a clear and steep velocity gradient across the cone, with average flux-weighted velocities from v < −200 km s⁻¹ at the edges to v > 100 km s⁻¹ at the centre along the outflow axis. Such a steep velocity gradient is commonly associated with a hollow conical geometry (Venturi et al. 2017; Venturi & Marconi 2020). The velocity dispersion of NGC 7582 and Circinus has deep minima along the galaxy plane, while it is larger along the outflow extension, with sporadic peaks of ∼350 km s⁻¹. NGC 4945 instead is characterised by a smoother and regular pattern, with the line width increasing to ∼400 km s⁻¹ along the outflow axis and then decreasing towards the edges down to ∼250 km s⁻¹, suggesting a more compact outflow geometry. We defined the systemic gas velocity in each source as the fitted stellar velocity, assuming co-rotation of gas and stars in the galaxy disc.

4.4. Outflow models

In this section, we present the results of our modelling applied to the three selected galaxies. Given the best set of parameters, we first created the corresponding model and provided spatially resolved estimates of the outflow energetics, and finally compared it with the volume-averaged results obtained with the standard literature method, described later in Sect. 4.5.1 (Cano-Díaz et al. 2012; Cresci et al. 2015b; Fiore et al. 2017; Harrison et al. 2018). From the data, we fixed the position angle (γ, measured counterclockwise from the north direction), the outflow conical aperture, and the centre of the outflow model. γ coincides with the projected approaching outflow axis; the outflow aperture and centre are instead identified by the ionised gas velocity and velocity dispersion maps. The outflowing cone aperture and centre are defined by the mentioned moment maps since these latter provide a better constraint on the projected shape and outflow starting point than the clumpy integrated emission. We then run the fit of the [OIII] observed emission line as described in Sect. 2, accounting for degeneracies and minimising the differences between modelled and observed spectra. To estimate the outflow velocity uncertainty, we fixed all the parameters except the outflow velocity, which was varied until measuring a variation from the minimum of the κ estimator of 10% (Eq. (6)). Figure 12 shows a comparison between the observed moment maps and those obtained with our best-fit models for NGC 4945 and NGC 7582; the observed data are masked to facilitate a comparison with the modelled maps. The same maps for Circinus are reported in first and third panels from the top in Fig. 13.

Fig. 12. Comparison between modelled moment maps and masked observations. From left to right, panels show the integrated flux, the LOS velocity and the velocity dispersion. Conical model for NGC 4945 (a) and NGC 7582 (b). Different models for Circinus are shown in Fig. 13.

Fig. 13. Comparison between observed data and different model configurations for Circinus galaxy. From left to right, the ionised [OIII]λ5007 emission line flux, the LOS velocity and velocity dispersion map. From top to bottom: masked observed data, best-fit model obtained with the parameters listed in Table 1, unweighted model, model with low maximum outflow radial velocity of 250 km s−1, and hollow conical model with inner opening angle of 27°. Each model map has the corresponding residual map on the right side. The maps are colour-coded according to the colour-bar scale of the observed data in the top panel.

The grainy emission and complex observed velocity fields of each outflow are extremely well reproduced by the model as a consequence of weighting clouds by the observed flux. Velocity dispersion maps highlight some discrepancies; for example in Circinus, where the modelled velocity dispersion does not properly reproduce the observed clumpy structure. There might be two possible explanations for this. First, the disc contribution may not be negligible, and therefore we should improve the model by considering the emission from background gas in the galaxy disc. Second, there could be intrinsic velocity dispersion in the outflow, which we have not considered. As reported in Table 1, we find outflows with axes close to the plane of the sky (i.e. β ∼ 90°) in all sources. This is an obvious consequence of the fact that we selected Seyfert-II galaxies with clear evidence for an extended ionisation cone, which implies that the LOS is close to perpendicular to the cone axis.

As shown by the best-fit models in Fig. 12 for NGC 4945 and NGC 7582 and in Fig. 13 for Circinus, outflow features are perfectly reproduced by a full conical geometry. However, as in many studies the best result is provided by a hollow conical geometry, we tested our model with hollow conical geometries (e.g., Crenshaw et al. 2000; Crenshaw & Kraemer 2000; Das et al. 2005; Bae & Woo 2016, for a statystical analysis).

Therefore, we ran a fitting procedure considering different inner opening angles to reproduce the data. As shown in Fig. 13 for Circinus, this configuration is not able to reproduce the data; the same conclusions apply to each galaxy in our sample. We show different model configurations compared to the observed moment maps for Circinus (top panel). For each model, the residual moment maps are also shown, which are obtained by subtracting the observed and model emission spaxel by spaxel. Both the hollow cone and the low outflow velocity configurations are clearly unable to reproduce the observed emission. Indeed, hollow cones miss clouds at small projected velocities, while the low-outflow-velocity model is unable to reproduce the observed line wings (see e.g., Fig. 5). As the cone fills up, the residuals improve, reaching a minimum when the cone is completely filled. We conclude that a full conical geometry –with slow changes in radial velocity around a mean constant value– is well suited to reproducing the observed outflow features in our sample.

4.5. Outflow energetics

In this section, we investigate the ionised mass and the energetics of the outflow using two different approaches: the first approach (hereafter, referred to as the ‘standard method’) follows the prescription of Cano-Díaz et al. (2012; also e.g., Marasco et al. 2020; Tozzi et al. 2021) and relies on several assumptions, because outflow physical properties are usually unknown. The second (hereafter referred to as the MOKA^3D based method) is based on our novel 3D geometry and kinematic modelling.

4.5.1. Standard method

We calculate the mass-outflow rate through the surface of a spherical sector with radius r defined by the cone aperture, using the outflow velocity inferred from the [OIII] emission line. If the outflow geometry and inclination with respect to the LOS are unknown, as done for example in Marasco et al. (2020) and Tozzi et al. (2021) we can define the outflow velocity as

(8)

where and are the maximum percentile velocities corresponding to 10% and 90% of the flux of the outflow component profile, respectively, and v_sys is the systemic velocity of the galaxy, set to 0 km s⁻¹. Equation (8) assumes that the intrinsic outflow velocity can be determined by the wing of the line profile, which provides the maximum observed velocity, possibly correcting for projection effects and dust absorption (e.g., Cresci et al. 2015b). The [OIII] line luminosity (L_[OIII], that we measure in each spaxel from the observations) can be written from a theoretical point of view as

(9)

where V is the volume occupied by the outflowing ionised gas, f the filling factor of the [OIII] emitting clouds in the outflow, n(O²⁺) and n_e are the volume densities of the O²⁺ ions and of electrons, respectively, and ϵ_[OIII] the [OIII]λ5007 emissivity, having a weak dependence on the temperature (∝T^0.1) at the typical temperature of the NLR (∼10⁴ K). As done in Cano-Díaz et al. (2012), here we assume that most of the oxygen within the outflowing ionised gas is in the O²⁺ form, and neglect the contribute to the mass from species heavier than helium. Finally, the outflow mass can be computed as

(10)

where L₄₄([OIII]) is the luminosity of the total [OIII]λ5007 emission line profile in units of 10⁴⁴ erg s⁻¹, ⟨n_e, 3⟩ is the average electron density in the ionised gas clouds in units of 10³ cm⁻³ (i.e. ⟨n_e, 3⟩=∫_Vn_e f dV/∫_V f dV) and [O/H] represents the oxygen abundance in solar units. is a ‘condensation factor’ (for more details on the derivation of the warm mass traced by the [OIII]λ5007 line in every pixel see Rupke et al. 2005; Genzel et al. 2011; Carniani et al. 2015; Cresci et al. 2017; Veilleux et al. 2020). We can assume C = 1 under the simplifying hypothesis that all ionised gas clouds in each resolution element (a MUSE spaxel in our case) have the same density. Also, under these assumptions, the mass of the outflowing ionised gas is independent from the filling factor of the emitting clouds. Finally, we are able to infer the average mass outflow rate in the volume, as follows:

(11)

where v₃ is the outflow velocity in units of 1000 km s⁻¹, and R_kpc is the conical outflow radius, in units of kpc. The outflow rate is independent of both the opening angle Ω of the outflow and the filling factor f of the emitting clouds (under the assumption of clouds with the same density). The kinetic energy (E_kin), kinetic luminosity (L_kin), and momentum rate (ṗ_out) of the ionised outflow, are then given by:

(12)

(13)

(14)

Equations (10)–(14) require knowledge of the different physical properties of the outflow, but only a few are usually measured, while others have to be assumed. The quantities usually assumed are the oxygen abundance, which is usually fixed to the solar abundance, and the electron density, which can be estimated from the [SII]λλ6717, 6731 doublet or needs to be fixed to typical values of AGN at a similar redshift to the sample. If the S/N for the [SII] doublet is high enough and the two lines are spectrally resolved, n_e can be directly estimated from the flux ratio for the lines in the doublet (e.g., Osterbrock 1989; Kewley et al. 2019), as follows:

(15)

where R is the flux ratio of the total emission line profile of the sulfur doublet f([SII]λ6717)/f([SII]λ6731). Assuming a constant electron density can have a huge impact on the outflow energetic (e.g., Veilleux et al. 2018; Davies & Baron 2020), but is nevertheless necessary in those cases where an estimate from the data is not possible.

The outflow energetic properties – calculated using Eqs. (10)–(14) – are reported in Table 4, except for the ionised mass outflow rate, which is discussed below. Their uncertainties were obtained with error propagation; for the electron density we assumed a systematic uncertainty of 50% (Tozzi et al. 2021). For NGC 7582 and NGC 4945, we determined the electron density from [SII] doublet with Eq. (15). For Circinus, instead we averaged the value provided by the spatially resolved study by Mingozzi et al. (2019), which was obtained from the high-velocity parts of the [SII] doublet lines. We find mass-outflow rates of 1.9 × 10⁻³ M_⊙ yr⁻¹, 0.2 × 10⁻² M_⊙ yr⁻¹, and 1.4 × 10⁻² M_⊙ yr⁻¹ for NGC 4945, Circinus, and NGC 7582, respectively. NGC 4945 is also characterised by a molecular outflow traced by ALMA CO J = 3–2 emission that is co-spatial to the ionised one, with an estimated mass-outflow rate of ∼20 M_⊙ yr⁻¹, as reported by Bolatto et al. (2021; Carniani et al., in prep.). The inferred Ṁ_out for Circinus, computed over a distance of ∼1 kpc from the AGN, is smaller by a factor of approximately 100 compared to the estimate of Fonseca-Faria (2021). These latter authors found outflow rates for the blue and red components of 0.12 M_⊙ yr⁻¹ and 0.1 M_⊙ yr⁻¹, respectively, by selecting the spaxels with detected high-ionisation outflow traced by [Fe vii]λ6089, and considering a maximum distance from the AGN of 700 pc. The Ṁ_out obtained for NGC 7582 is consistent with the high-ionisation counterpart of 0.7 × 10⁻² M_⊙ yr⁻¹ inferred by Davies & Baron (2020) using VLT/X-shooter spectra and covering the inner 300 pc. Both the outflow extension and maximum velocity of 364 km s⁻¹ are much lower compared to our work, in which we estimated the outflow rate considering an outflow extension of ∼3.2 kpc and a radial velocity of ∼630 km s⁻¹.

4.5.2. Wind energetic uncertainties

Literature estimates of the wind parameters rely on different assumptions for each gas phase. The electron density in the outflowing gas, the temperature, and geometry are among the parameters that most affect the uncertainties of the outflow energetics. As explained in Sect. 4.5.1, when it comes to estimating the ionised wind mass and energetics, many basic assumptions can affect the final parameters, leading to systematic uncertainties of a factor of about 10. Therefore, even when adopting the standard recipe outlined in Sect. 4.5.1, it is not unexpected to observe up to three orders of magnitude of scatter for outflow energetics in different studies (e.g., Fiore et al. 2017; Cicone et al. 2018; Harrison et al. 2018). In this context, our MOKA^3D method represents an innovative approach to determining the outflow physical properties with great accuracy, because it is not based on strong observational or physical assumptions.

4.5.3. MOKA^3D based method

Employing the tomographic outflow reconstruction, we have a 3D distribution of ionised gas clouds within the outflow. This allowed us to provide a spatially resolved estimate of the outflow energetics. To compute the amount of ionised mass in each voxel, we used Eq. (10), converting the flux density emitted by each cloud in the same bin to ionised mass in solar masses. Assuming the flux density within each spaxel to be constant over time and the outflow to subtend a solid angle dΩ, we use the continuity equation in spherical coordinates and express the ionised mass outflow rate via the following equation:

(16)

where ρ is the mass density in each 3D volume element and v is the outflow velocity, which is assumed to be constant within each spatial element. The mass density in each voxel is ρ = dM_ion/dΩr²dr, with dr = 0.2″ being the MUSE spatial resolution and dM_ion the ionised mass in the same spaxel. Therefore, Eq. (16) becomes

(17)

where Ṁ_out is the average mass-outflow rate within each spatial element of width dr at a certain distance from the AGN. To estimate the mass-outflow rate profile, we computed the amount of ionised gas (dM_ion) in a shell of fixed width (dr) and assumed the radial velocity to be constant within the shells. Two-dimensional maps of the mass-outflow rate and their respective radial profiles, for NGC 4945, Circinus, and NGC 7582 are shown in Fig. 14 in the left and right panels, respectively. To calculate the uncertainties on Ṁ_out in each shell, we propagated the errors, considering the inferred error on v from our modelling (see Sect. 4.4) and assuming an uncertainty of 0.6″ on dr.

Fig. 14. Mass outflow rate distribution of the ionised gas derived from [NII]λ6584 for NGC 4945 and from [OIII]λ5007 for Circinus and NGC 7582. Left panels: 2D maps of the mass-outflow rate distribution. Right panels: radial profile of the mass-outflow rate from the cone vertex. The volume average mass-outflow rate obtained with the MOKA3D-based method (dashed green line) is compared to the standard method result (dotted blue line). The standard method estimate of the volume-average Ṁout is underestimated by a factor of two in NGC 4945 and Circinus. From top to bottom: NGC 4945, Circinus, and NGC 7582.

We find three different radial profiles, which can be linked to the past AGN accretion histories and variation of AGN luminosity, possibly indicating large time variations in ionising luminosity that result in the peaks of the outflow rate (Keel et al. 2012, 2015, 2017; Gagne et al. 2014; Finlez et al. 2022). However, this could also be an ionisation effect not taken into account when assuming a constant luminosity-mass conversion factor, as we have done. A proper determination of the ionised gas mass from the emission-line luminosity requires the use of photo-ionisation models, which we will address in future work. This can also be observed by inspecting the outflow clumpiness distribution. In NGC 4945 and Circinus, we observe two shells with increased flux density. In NGC 7582 instead, we observe one peak of emission slowly decreasing with the distance from the source, up to the maximum model extension. From the mass-outflow rate profile, we can estimate the dynamical timescale of the ionised outflow t_dyn = d_e/v_e, where d_e and v_e are the maximum outflow extension and intrinsic velocity inferred with MOKA^3D, respectively. We obtained 1.4 × 10⁶ yr, 1.8 × 10⁶ yr and 4.9 × 10⁶ yr, for NGC 4945, Circinus, and NGC 7582 respectively, which are consistent with theoretical predictions and simulations (Nardini & Zubovas 2018; Richings & Faucher-Giguère 2018).

4.5.4. Methods comparison

In addition to the spatially resolved mass-outflow rate estimate, MOKA^3D can provide a volume-averaged estimate of Ṁ_out. This is done using Eq. (17) by considering the total amount of ionised gas mass inferred by the emission of each model cloud, the maximum model distance from the AGN, and the intrinsic outflow velocity. As explained in Sect. 2.4, each cloud has an assigned weight which is derived from the observed emission in the corresponding voxel. Therefore, we convert each cloud weight to [OIII] luminosity and then derive the ionised gas mass. Table 5 shows the parameters and the volume-average mass-outflow rate estimated with the standard method (Sect. 4.5.1) and our MOKA^3D model (Sect. 4.5.3). The estimates provided by the two methods are compatible within the errors reported in Table 5. Moreover, the volume-average mass-outflow rates have the same order of magnitude as the radial profile, as shown in the right panels in Fig. 14. The standard method underestimates the average Ṁ_out by a factor of 2 compared to our method based on MOKA^3D. This is due to the fact that both the velocity (see Eq. (8)) and the projected maximum radius extension adopted by the standard method are smaller than the values provided by our model. Both methods are affected by systematic uncertainties of the order of 50% as a result of the assumed uncertainties on the electron density. As discussed above, to refine the energetics estimate, we plan to combine our MOKA^3D model with a photo-ionisation model in order to properly constrain the outflowing gas mass and further reduce the uncertainties.

5. Conclusions

We developed a new 3D kinematical model that allows us to constrain the kinematics and geometry of clumpy AGN galactic winds with unprecedented accuracy, reproducing the emission line features observed in all spaxels. We applied the model to three nearby Seyfert-II galaxies selected from the MAGNUM survey (Cresci et al. 2015b; Venturi et al. 2018, 2021; Mingozzi et al. 2019), featuring a clear (bi)conical ionised outflow extending over kpc scales, and showed that the observed complexity in the kinematical maps is well reproduced by a simple radial outflow model with constant velocity and a clumpy ionised gas cloud distribution. The main features and results of our model are summarised below:

1. In Sect. 2, we present a multi-cloud model that takes into account the intrinsic clumpy nature of ionised outflows. We weight each of our model clouds based on the observed line emission and velocity, and constrain the input parameters through a fit by comparing the model and observed line profile spaxel by spaxel.

2. The model takes into account the spatial and spectral resolution of the data by creating a model cube with identical 3D extensions and spatial and spectral resolutions of observed data cubes. MOKA^3D also accounts for the convolution of the surface brightness of model clouds with the PSF measured from observed data.

3. One of the main achievements of this modelling technique is the disentanglement of the observed features from projection effects, obtaining the real 3D distribution of the gas clouds in a tomographic way.

4. In this work, the free outflow kinematical and geometrical parameters that the modelling retrieves are the outflow inclination with respect to the plane of the sky, the gas cloud radial velocity, and the systemic velocity with respect to the host galaxy (see Sect. 4.4).

5. In Sect. 3, we tested the capabilities of our new method by applying it to four simulated cubes with moment maps similar to the ones observed in the MAGNUM sample. The model manages to reproduce the observed features and provides the correct outflow parameters, with an accuracy > 95%, fitting up to four parameters.

6. In Sect. 4.5.3, we show how we can measure the outflow physical properties (i.e. mass outflow rate, kinetic luminosity, kinetic energy, and momentum rate) with great accuracy and without strong assumptions.

With our model we provide a new and reliable method to constrain the wind energetics and the powering mechanism in AGN, and a different approach to inferring the impact of outflows on the ISM and galaxy evolution, expanding our knowledge of AGN feedback mechanisms. In Cresci et al. (2023), we present the first application of MOKA^3D to a high-z QSO observed with JWST/NIRSpec and characterised by a collimated high-velocity outflow piercing an expanding bubble of ionised gas³. In future work, we will apply our model to the other galaxies in the MAGNUM survey, as well as to other high-redshift quasars (e.g., those in Carniani et al. 2015). We will also consider more complex models where a disc or other kinematical components contribute to the observed data cubes.

Acknowledgments

G.C., A.M., G.T., F.M., F.B. and G.V. acknowledge the support of the INAF Large Grant 2022 “The metal circle: a new sharp view of the baryon cycle up to Cosmic Dawn with the latest generation IFU facilities”. G.C. and A.M. acknowledge support from PRIN MIUR project “Black Hole winds and the Baryon Life Cycle of Galaxies: the stone-guest at the galaxy evolution supper”, contract # 2017PH3WAT. E.D.T. was supported by the European Research Council (ERC) under grant agreement no. 101040751. S.C. acknowledges funding from the European Union (ERC, WINGS,101040227). G.V. acknowledges support from ANID program FONDECYT Postdoctorado 3200802.

References

All Tables

All Figures

Fig. 1. Flowchart of our MOKA3D kinematic model. The method consists in assigning the model parameters and assuming a coordinate system and a velocity field. Each cloud is assigned the coordinates in the source reference frame (green), which is later transformed to the observer reference frame (blue), through Euler rotation matrices. If necessary, the method requires the convolution with the PSF and/or the spectral LSF. Then, the model is binned in the space sampled by the observed data cube (red). The model emission line profile percentiles at 1% and 99% are compared to the observed ones and, if they are comparable, each particle is weighted according to the observed emission in the corresponding bin. The model is evaluated by means of the κ estimator (Eq. (6)). The free parameters are varied and the loop repeats, until finding the suitable parameter configuration that minimises κ and reproduces the observed spectrum in each spaxel.

In the text

	Fig. 2. Source (blue) and observer (red) reference frame. Euler angles α, β, and γ are used to transform from the source frame to the observer’s frame and vice versa. The purple dotted line represents the line of nodes, N. The observer z axis coincides with the LOS, and (x, y) is the plane of the sky.
In the text

Fig. 3. Coordinate transformation and binning procedure of the model: from the simple conical geometry in the source reference frame, uniformly populated with clouds, to the 3D model cube in the observer’s reference frame. The model is then binned to assign to each cloud a position in the (x, y, vz) space, as highlighted for the purple cloud. The model cube has the same spatial and spectral resolution as the observed data cube in order to weight the clouds according to the observed flux in the corresponding spaxel. Dashed grey lines show two adjacent voxels, corresponding to different spatial and/or spectral channels.

In the text

Fig. 4. Example of unweighted models simulating an outflow with the spatial and spectral resolution of a galaxy observed with MUSE at z = 0.001. From left to right: Line-velocity integrated map, LOS velocity map, and velocity dispersion map. Top panels: full biconical model with P.A. = 40°, inclination with respect to the plane of the sky of β = 85°, a constant radial velocity of 1200 km s−1, and a cone semi-aperture of θOUT = 30°. Bottom panels: single hollow conical model with P.A. = 90°, β = 65°, and inner and outer semi-aperture θIN = 25° ,θOUT = 35°. Both models have been convolved with a spatial PSF = 0.3″ and no spectral convolution.

In the text

Fig. 5. Comparison between the observed [NII]λ6584 emission line profile (black) and the model emission (red) of NGC 4945, extracted from a spaxel with a prominent blue-shifted wing. Top and bottom panels show the unweighted and weighted emission, respectively. Left panels: the model emission is obtained with the best fit parameters reported in Table 1. Right panels: as a demonstrative example, we show the model line profile as obtained by assuming an unsuitable inclination and opening angle. Such a configuration is clearly not able to to reproduce the high-velocity, blue-shifted emission, thus causing a sharp cut off. The absence of clouds in all velocity channels where emission is observed, prevent the model to reproduce the observed outflow properties.

In the text

Fig. 6. Comparison of observed data and two degenerate models of NGC 4945. Moment maps of observed data (top panels) and of two models with high radial outflow velocity and incorrect inclination with respect to the plane of the sky, incorrect outer opening angle, and incorrect systemic velocity (middle and bottom panels) with respect to the best model shown in Fig. 12. See Sect. 2.6 for details of the parameters of the degenerate models and how to overcome them.

In the text

Fig. 7. Schematic explanation of the effect of an excessive radial outflow velocity with fixed model aperture and inclination. The solid black lines define the outflowing cone aperture, the dotted black lines show the LOS directed towards the observer, and the grey dotted lines show the region containing the clouds with projected velocity within the v1, v99 observed boundaries. The grey clouds have LOS velocity above v1, v99, i.e. their weights are zero. Black arrows show the decomposition of the outflow velocity. The top panel shows the view from the top of the model, while the bottom panel shows the flattening effect as seen from the plane of the sky.

In the text

Fig. 8. Moment maps and integrated emission comparison between three simulated and best-fit models. From top to bottom: Model 1, Model 2, and Model 3, whose parameters are listed in Table 2. Left panel: comparison of simulated and best-fit moment maps. From left to right: integrated emission, LOS velocity, and velocity dispersion maps. Right panel: comparison of simulated data (dotted blue) and best-fit integrated profile (red). Residual emission is shown in green.

In the text

	Fig. 9. Same as Fig. 8 but for Model 4. Left panels: moment maps comparison of simulated data (top) and best-fit model (bottom) obtained with four free parameters. From left to right: integrated emission, LOS velocity, and velocity dispersion maps. Right panel: integrated emission line comparison between simulated (dotted blue) and best-fit model (red), with residuals in green.
In the text

Fig. 10. Example of outflow- and disc-selected clumps in the Circinus galaxy. Top panel: 2D [OIII]λ5007 emission map of the Circinus galaxy. Orange and black circles show the average regions where the disc- and outflow-integrated emission are extracted, respectively. Bottom panel: comparison of outflow and disc-integrated spectra (×102.5) extracted from the corresponding regions shown in the left panel.

In the text

Fig. 11. Moment maps of the sample. From top to bottom: NGC 4945, Circinus, and NGC 7582. From left to right: Integrated surface brightness map in units of 10−20 erg cm−2 s−1 arcsec−2, intensity-weighted velocity and velocity dispersion. The first row shows the [NII]λ6584 emission in NGC 4945 and the last two rows refers to the [OIII]λ5007 emission. North is up. The maps are 1 spaxel−σ spatially re-smoothed, so as to get a better visual output, and pixels with S/N ≥ 3 are masked. The maps are obtained from the spectroscopic analysis of the modelled emission line from MUSE data (e.g., Marasco et al. 2020).

In the text

	Fig. 12. Comparison between modelled moment maps and masked observations. From left to right, panels show the integrated flux, the LOS velocity and the velocity dispersion. Conical model for NGC 4945 (a) and NGC 7582 (b). Different models for Circinus are shown in Fig. 13.
In the text

Fig. 13. Comparison between observed data and different model configurations for Circinus galaxy. From left to right, the ionised [OIII]λ5007 emission line flux, the LOS velocity and velocity dispersion map. From top to bottom: masked observed data, best-fit model obtained with the parameters listed in Table 1, unweighted model, model with low maximum outflow radial velocity of 250 km s−1, and hollow conical model with inner opening angle of 27°. Each model map has the corresponding residual map on the right side. The maps are colour-coded according to the colour-bar scale of the observed data in the top panel.

In the text

Fig. 14. Mass outflow rate distribution of the ionised gas derived from [NII]λ6584 for NGC 4945 and from [OIII]λ5007 for Circinus and NGC 7582. Left panels: 2D maps of the mass-outflow rate distribution. Right panels: radial profile of the mass-outflow rate from the cone vertex. The volume average mass-outflow rate obtained with the MOKA3D-based method (dashed green line) is compared to the standard method result (dotted blue line). The standard method estimate of the volume-average Ṁout is underestimated by a factor of two in NGC 4945 and Circinus. From top to bottom: NGC 4945, Circinus, and NGC 7582.

In the text