A disc wind model for blueshifts in quasar broad emission lines

ABSTRACT

Blueshifts – or, more accurately, blue asymmetries – in broad emission lines such as C iv λ1550 are common in luminous quasars and correlate with fundamental properties such as Eddington ratio and broad absorption line (BAL) characteristics. However, the formation of these blueshifts is still not understood, and neither is their physical connection to the BAL phenomenon or accretion disc. In this work, we present Monte Carlo radiative transfer and photoionization simulations using parametrized biconical disc-wind models. We take advantage of the azimuthal symmetry of a quasar and show that we can reproduce C iv blueshifts provided that (i) the disc-mid-plane is optically thick out to radii beyond the line formation region, so that the receding wind bicone is obscured; and (ii) the system is viewed from relatively low (that is, more face-on) inclinations (≲40°). We show that C iv emission-line blueshifts and BALs can form in the same wind structure. The velocity profile of the wind has a significant impact on the location of the line formation region and the resulting line profile, suggesting that the shape of the emission lines can be used as a probe of wind-driving physics. While we are successful at producing blueshifts/blue asymmetries in outflows, we struggle to match the detailed shape or skew of the observed emission-line profiles. In addition, our models produce redshifted emission-line asymmetries for certain viewing angles. We discuss our work in the context of the C iv λ1550 emission blueshift versus equivalent-width space and explore the implications for quasar disc wind physics.

1 INTRODUCTION

Luminous active galactic nuclei (AGNs) and quasars produce outflows in the form of collimated, relativistic jets and mass-loaded winds. Both of these classes of outflows are intimately connected to the accretion disc and have profound implications for our understanding of black holes, particularly their growth and evolution. In addition, outflows are invoked as possible feedback agents (see Fabian 2012; Morganti 2017; Harrison et al. 2018, for reviews), providing a mechanism through which the vast accretion energy released close to the black hole can couple to the surroundings.

AGN winds – our focus here – are capable of expelling gas and/or quenching star formation in their host galaxies (e.g. Silk & Rees 1998; King 2003; Di Matteo, Springel & Hernquist 2005; Feruglio et al. 2010; Hopkins & Elvis 2010; Costa, Sijacki & Haehnelt 2014), and are critical ingredients in cosmological simulations (e.g. Beckmann et al. 2017; Chisari et al. 2018). In addition, the impact of these winds can feasibly account for observed scaling relations such as the M–σ relation (Silk & Rees 1998; King 2003; Di Matteo et al. 2005). Despite the clear importance and apparent ubiquity of outflowing material, the driving mechanism of the winds and the precise relationship between the outflows and the accretion process both remain poorly understood.

Approximately 20−40 per cent of quasars show blueshifted broad absorption lines (BALs) in their rest-frame ultraviolet (UV) spectra (Weymann et al. 1991; Hewett & Foltz 2003; Knigge et al. 2008; Allen et al. 2011; Dai, Shankar & Sivakoff 2012). These BAL quasars provide unambiguous evidence of an absorbing outflow along the line of sight. The general paradigm is that BALs are produced by an accretion disc wind that intersects the observer’s sightline. A fairly common geometry is an equatorial wind (e.g. Murray et al. 1995), roughly consistent with simulations of line-driven winds (Proga & Kallman 2004). However, this equatorial disc-wind geometry is hard to reconcile with the apparent similarity in spectral properties between BAL and non-BAL quasars (e.g. Weymann et al. 1991; Matthews, Knigge & Long 2017; Yong et al. 2018; Rankine et al. 2020). Furthermore, an equatorial viewing angle would lead BALs to be under-represented in quasar samples due to suppression of the continuum flux by an obscuring torus, foreshortening of the accretion disc and/or limb darkening (Krolik & Voit 1998; Matthews et al. 2017). In addition to questions of geometry, measuring distances to BAL formation regions throws up further puzzles – how can apparent distances of ∼kpc scales (in many cases) inferred from density sensitive lines (e.g. Borguet et al. 2012; Arav et al. 2018; Xu et al. 2018; Leighly et al. 2022) be reconciled with the much smaller launching radii expected from line-driven disc wind simulations (e.g. Proga & Kallman 2004)? The former inferred BAL distances are much larger than the scale of the quasar accretion disc, which has a self-gravitation radius of ∼0.1 pc. All in all, although BALs must be formed in quasar outflows of some form, whether these outflows are actual ‘disc winds’ in the conventional sense is not at all clear. Indeed, perhaps BAL quasars are not a single, monolothic class of objects, and the physics driving their observational appearance could be as heterogeneous as the spectral signatures themselves.

Although blueshifted BALs are the ‘smoking gun’ signature of an absorbing outflow, a disc wind can also affect the spectrum in other ways. For example, various authors have proposed a common origin for the formation of broad emission lines and BALs within a biconical wind structure (e.g. Emmering, Blandford & Shlosman 1992; de Kool & Begelman 1995; Murray et al. 1995; Elvis 2000). Emission-line asymmetries can also be produced in outflows, but these are more challenging to interpret than BALs. While we know that the emission from the blue wing comes from material moving towards us, this is subtly different from knowing it is formed in an outflow. Assuming symmetry about the disc axis, an outflow also has a receding emitting volume. Whereas in a pure scattering atmosphere with spherical symmetry a classical P-Cygni profile is produced, for the more complex quasar system the observed profile depends on a number of factors; these factors include the covering factor of the flow, the emissivity profile of the line, the opacity of the mid-plane and the interplay between inclination and wind launching angle. Indeed, the early papers describing the discovery of C iv λ1550 emission line blueshifts discussed formation in a quasar outflow (Gaskell 1982; Wilkes 1984), but alternative models in which the BLR is inflowing have also been proposed (Gaskell & Goosmann 2013, 2016).

Extending the ‘Eigenvector I’ paradigm (Boroson & Green 1992; Sulentic et al. 2000; Marziani et al. 2001; Sulentic et al. 2002; Marziani et al. 2003) to higher redshift, Sulentic et al. (2007) and Richards et al. (2011, hereafter R11) identified the fundamental nature of the C iv equivalent width (EW) versus blueshift parameter space (hereafter C iv emission space). R11 demonstrated the existence of fundamental trends across this parameter space in terms of emission-line EWs, Eddington ratio, and radio loudness. Rankine et al. (2020, hereafter R20) expanded on this by using mean-field independent component analysis to reconstruct BAL and non-BAL quasar spectra, showing that the BAL quasars appear to be drawn from the same parent population as non-BAL quasars. The results presented by R11 and R20 have important implications for our understanding of quasar outflows; they show that the strongest, broadest BAL troughs are formed in quasars with high C iv blueshifts. The results can be interpreted within the framework of line-driven winds (Proga, Stone & Kallman 2000; Proga & Kallman 2004; Proga 2005), with high blueshift quasars corresponding to high Eddington fraction quasars with soft SEDs that provide the optimal conditions for line-driving (Rivera et al. 2022; Temple et al. 2023). However, it is important to remember that C iv blueshift versus EW space is likely a projection of a multidimensional space, and there are many questions that arise – for example, why are BALs only seen in most, but not all, of this parameter space? What is the role of inclination? How are the C iv blueshifts formed in the first place?

Radiative transfer calculations with varying degrees of complexity have been successful in simulating emission lines at observable levels from quasar disc wind models (Murray et al. 1995; Matthews et al. 2016; Waters et al. 2016; Matthews et al. 2020; Naddaf & Czerny 2022). A few authors have also specifically discussed synthetic emission-line asymmetries and blueshifts in more detail. For example, Chajet & Hall (2013, 2017) used the magneto-hydrodynamic (MHD) wind model of Emmering et al. (1992), a slightly modified version of the classic Blandford & Payne (1982) self-similar solution, to study line widths and shapes as a function of the wind geometry. Yong et al. (2017) also studied line asymmetries using the same kinematic prescription from Shlosman & Vitello (1993) adopted here, including a simplified treatment of radiative transfer. They found they could produce blueshifts in their models, and that higher blueshifts were produced for polar winds and polar viewing angles due to velocity projection effects and the obscuring effect of the disc mid-plane. Although the observed blue asymmetries are most common in the C iv λ1550 line, they can be seen in other emission lines. For example, GRAVITY Collaboration et al. (2020) reported an asymmetric Paschen-α line in IRAS 09149-6206, also noting the importance of an opaque mid-plane in their modelling.

An interesting point of astrophysical comparison is protoplanetary systems, in which relatively low-temperature discs of gas and dust form around low-mass stars early in their lives (Williams & Cieza 2011). In these objects, blueshifted emission lines are observed in forbidden lines such as [Ne ii] 12.81 |$\mu$|m and O i λ6300. These blueshifted profiles are attributed to line formation in a disc wind, which, in T Tauri systems, is thought to be photoevaporative or thermally driven (Pascucci et al. 2020), although there may also be an inner MHD wind. The temperatures and kinematics of these winds are rather different from the putative broad-line region (BLR) or BAL winds, with typical velocities of ∼10 km s⁻¹ and electron temperatures of ∼3000 K. (Banzatti et al. 2019). The thermal winds are produced when the sound speed exceeds the local escape speed, and – though the heating and cooling physics is different – the winds bear similarity to the Compton-heated thermal winds thought to occur in X-ray binaries (Begelman, McKee & Shields 1983; Done, Tomaru & Takahashi 2018; Higginbottom et al. 2018) and possibly AGN (Mizumoto et al. 2019; Ganguly et al. 2021). A number of simulations of thermal/photoevaporative winds in protoplanetary discs have been successful in producing blueshifted profiles comparable to those observed (e.g. Sellek, Clarke & Booth 2021). In these simulations, the mid-plane of dust and cold gas is responsible for obscuring the receding wind bicone; this scenario is very similar to the one we will discuss for blueshift formation in quasar winds.

In this paper, we expand on previous efforts to model quasar emission lines (and their asymmetries) by conducting MCRT and photoionization simulations in an azimuthally symmetric biconical disc wind geometry. Our approach allows us to include a more physical and self-consistent treatment of the line formation mechanism compared to previous studies. We are also able to calculate the emissivity of the wind plasma, which can be used to infer information about the power and mass-loss rates of a disc wind associated with a C iv blueshift and study where the emission lines are actually forming in the wind. Our aim is not to model a population of quasars, or to reproduce the full diversity of quasar spectra. Instead, we will investigate whether C iv blueshifts can be produced for reasonable model assumptions within this assumed geometry, and whether they can explain the ‘windy’ component that is apparently needed to explain the diversity of quasar spectra. Our paper is structured as follows. We first introduce the observational properties of the C iv emission line and our numerical method in Section 2. Following an heuristic demonstration of our approach in Section 3, we present a series of simulated synthetic spectra in Section 4 and explore trends with the physical parameters and geometry of the wind model. In Section 5, we discuss our results in the context of quasar feedback, disc wind theory and line asymmetries in a variety of astrophysical settings, before concluding in Section 6. Throughout, we follow convention by using the shorthand C iv λ1550 to refer to the C iv doublet (corresponding to the resonant 1s²2p→1s²2s transitions) with rest frame vacuum wavelengths 1548.19 and 1550.77 Å (according to the NIST Atomic Spectra Database; Kramida et al. 2022).

2 OBSERVATIONAL DATA AND NUMERICAL METHODS

2.1 Observed Spectra

Before introducing our numerical methods and discussing the physics of line formation in disc winds, we first review what is known observationally about emission-line asymmetries. As stated above, C iv blueshifts were discovered some four decades ago (Gaskell 1982). We now know that, for typical quasars at z ∼ 1.5−3.5, the C iv λ1550 emission-line profiles form a continuum from strong and symmetric to weak, highly asymmetric and blueshifted.

The results of the simulations in the paper are presented primarily as C iv λ1550 emission-line profiles. As motivation and orientation for the investigation we use the spectroscopic sample from the Sloan Digital Sky Survey (SDSS) employed in R20. In particular, we will use this sample as the basis for presenting composite spectra and plotting the distribution of C iv λ1550 EWs and blueshifts.

To quantify the asymmetry in the emission line profiles, we follow the recipe of Coatman et al. (2016) to measure the EW and define blueshift of a line as

Here, λ_half is the wavelength that bisects the integrated line flux, and λ₀ is the line centre wavelength (1549.5 Å for the C iv UV resonance line doublet). This definition is positive if the line is blueshifted, opposite to the usual radial velocity definition. We also define the ‘skew’ of the line profile as the relative offset between the velocity of the emission-line peak and the blueshift, i.e.

where λ_peak is the wavelength of the peak line flux. We will use these two quantities to quantify of the asymmetry in our simulations and compare to the observed data, while using the skew as a measure of the line shape. This definition of skewness leads to a negative value when the line peak is redwards of λ_half.

Fig. 1 (right-panel) shows the C iv emission space – that is, the observed distribution of blueshift and EW for the C iv emission line – in the quasar sample from R20. For five locations in this parameter space, chosen so as to span a range of C iv emission-line properties, the left-hand panel shows high signal-to-noise ratio composite spectra (each generated using a minimum of 500 quasar spectra). The full set of 33 composite spectra, from which these five are drawn, are described by Stepney et al. (2023). The significant range of line profiles for C iv λ1550 and the C iii] λ1909 emission-line complex¹ is striking. The trend originally reported by R11 can be clearly seen: the high-blueshift quasars have weak He ii λ1640 emission and weaker C iii] λ1909 compared to the low-blueshift ‘peaky’ spectra. Additionally, the highly blueshifted, lowest EW C iv line in the composite spectra in Fig. 1 has a negative skew. In fact, negative skews when blueshifts are present in the C iv λ1550 are a general feature of SDSS quasar spectra (e.g. R11, R20; Stepney et al. 2023; see also Section 4.2.1).

2.2 MCRT and photoionization

We use a well-tested Monte Carlo radiative transfer (MCRT) and photoionization code, confusingly called python, to simulate rest-frame UV spectra from a biconical disc wind geometry. The code was originally introduced by Long & Knigge (2002), but has since been updated and upgraded (Sim 2005; Higginbottom et al. 2013; Matthews 2016; Matthews et al. 2016). Radiative transfer is treated using the Monte Carlo method, with line transfer treated in the Sobolev limit (see Noebauer & Sim 2019, for a review). During a simulation, the code follows energy packets of discretized radiation (hereafter ‘r-packets’, following the notation of Lucy 2002), as they make their way through a user-defined geometry or imported grid. The frequency distribution of these packets and their overall energy normalization are set by the input SED (Section 2.2.3). The code first runs a series of ‘ionization cycles,’ during which Monte Carlo estimators are used to iteratively calculate the temperature, ionization state, and emissivity of the wind. Once a satisfactorily converged solution is obtained (see Section 4), the temperature and ionization structure are held fixed. The synthetic spectrum is then created over a limited wavelength range (‘spectral cycles’) for a number of specified observer viewing angles.

2.2.1 Recent code improvements

Our MCRT code has been improved in a number of ways since it was last used to simulate quasar spectra. These improvements make only small quantitative differences to the spectrum observed from a quasar model. The first development since our work in Matthews et al. (2020, hereafter M20) is that we have improved our treatment of frame transformations and Doppler shifts so as to fully account for special relativistic effects. In particular, we now explicitly conserve energy in the co-moving frame and include the Lorentz factor in, for example, the calculation of the number density in a given volume. We also include all relevant special relativistic effects (such as relativistic aberration) when calculating photon frequencies.

An additional improvement has been made to the macro-atom mode of operation, in which transition probabilities are sampled to determine how energy flows into and out of the excitation, ionization, and kinetic energy pools. Previously, this was achieved by Monte Carlo sampling of the transition probabilities as originally proposed by Lucy (2002, 2003). Now, we are able to dramatically speed up this section of the code by using standard Markov chain techniques (Ross 2014), as applied to the MCRT macro-atom scheme by Ergon et al. (2018).

2.2.2 Wind model and clumping

As in our previous works, we use a flexible parametrized model for a biconical wind, originally introduced by Shlosman & Vitello (1993) as a model for disc winds in accreting white dwarf system. In this prescription, the wind rises from an accretion disc with straight poloidal streamlines between minimum and maximum radii, r_min and r_max, with a specified velocity law and range of opening angles. The velocity at a given distance l along a streamline is given by

and is controlled by four parameters: the acceleration length, R_v, the acceleration exponent, α, and the initial and terminal velocities, v₀ and v_∞. The rotational velocity, v_ϕ, is initially set to the Keplerian velocity at the streamline base, and specific angular momentum is then conserved along streamlines. The angles that the poloidal streamlines make with the vertical z-axis are spaced evenly between the inner and outer opening angles of the wind, θ_min and θ_max, respectively. A mass-loss rate, |$\dot{M}_w$|⁠, is imposed and the local density follows from mass conservation. The model is azimuthally symmetric around the polar axis, and has reflective symmetry around the disc mid-plane. We conduct our simulations on a cylindrical (x, z) grid with the wind discretized into logarithmically spaced cells. The number of cells is chosen to adequately resolve the dense wind base, and ensure cells do not have excessively large velocity changes or continuum optical depth.

As in our previous work, we use a simple treatment of clumping, inspired by modelling of stellar winds, known as ‘microclumping’ (e.g. Hamann, Koesterke & Wessolowski 1995; Matthews et al. 2016). Clumping is needed in order to prevent overionization of the flow; without it, emission lines are only formed if either the X-ray source is weaker than expected for a quasar of this luminosity (Higginbottom et al. 2013; Matthews et al. 2016), or the mass-loss rate through the wind is significantly higher than the accretion rate through the disc. In fact, we have to make a few assumptions in order to get the wind to emit strongly at all, which we discuss further in Sections 4.2.2 and 5.5.

2.2.3 Radiation sources and sinks

The r-packets in our MCRT modelling are all launched from a central, isotropic source, with a frequency distribution that matches the SED used by M20. This SED shape is a combination of a multitemperature, blackbody accretion disc with a standard Shakura & Sunyaev (1973) temperature profile, and a high-energy power law with spectral index 0.9 (that is, F_ν ∝ ν^−0.9). We explore the sensitivity of our results to the adopted SED and its anisotropy in Section 4.5, including a simulation in which the disc component of the radiation field is anisotropic and samples the emissivity profile of a multitemperature accretion disc. The r-packets are propagated through the simulation grid until they escape the system or strike an opaque surface; if the latter occurs, the r-packet is discarded and any heating effect is not accounted for in this work. The inner 6r_g (where r_g is the gravitational radius) is always a purely absorbing opaque surface in our simulations. This choice corresponds to the innermost stable circular orbit of a non-spinning black hole. The mid-plane is either set to be transparent or opaque, with a parameter R_mid that sets the horizontal distance at which this mid-plane becomes transparent again. An opaque mid-plane is crucial for forming the blueshifts we observe.

3 HEURISTIC DEMONSTRATION

We begin by using a single MCRT simulation of a biconical disc wind to demonstrate the basic principle underlying the blueshift formation in our models. We use a model comparable to Model A from M20, with full model parameters given in Table 1. Broadly speaking, the model is designed to mimic a typical quasar in terms of its mass and luminosity. The wind has a mass-loss rate equal to the accretion rate and is clumpy, with a filling factor of f_V = 0.01. The wind is launched between 300 and 600r_g, and has an inner wind opening angle of θ_min = 45°, a covering factor of Ω/4|$\pi$| = 0.2, and an acceleration exponent of α = 1. This model is denoted model n in subsequent discussions and figures.

Fig. 2 shows a representation of how blueshifted asymmetries can form in this axisymmetric system. The figure is intended as heuristic or illustrative, but it is not a schematic diagram in the sense that the plotted colourmaps and synthetic spectra are obtained from the full, self-consistent MCRT photoionization calculation. The wind geometry is shown, colour-coded by projected velocity for an observer at θ_i = 15°. We discuss the inclination dependence further in Section 4.1. The C iv λ1550 line profile is also plotted, and shows the results with and without an obscuring mid-plane disc (OMD). A clear blueshifted asymmetry can be seen in the opaque mid-plane case, as the receding wind cone is obscured. There is a slight redshifted asymmetry in the case with a transparent mid-plane due to slight preferentially absorption of the blue wing by the wind itself.

Fig. 2 provides a heuristic, but quantitative, demonstration of the basic principle of obscuring the receding part of the wind to break the symmetry of the system (and of the line profile). This general mechanism for forming line profile asymmetries – obscuration of the receding part of a quasar outflow, possibly by an extended disc-like structure – was discussed even in the earliest studies of C iv blueshifts (Gaskell 1982; Wilkes 1984). Since then, this scenario has been explored in more detail using radiative transfer, albeit with simpler methods than employed here (Chajet & Hall 2013, 2017; Yong et al. 2017; GRAVITY Collaboration et al. 2020), and similar proposals have been made for the line asymmetries in protoplanetary discs (Sellek et al. 2021). A related model was discussed by Richards et al. (2002), in which the mid-plane was only opaque at high inclinations due to limb opacity/path length effects. Similarly, blueshifts in AGN narrow lines such as O iii λλ4960,5008 have been attributed to, and modelled with, biconical winds on larger (∼100s of pc) narrow-line region scales with circumnuclear or mid-plane dust obscuring the receding bicone (Crenshaw et al. 2010; Bae & Woo 2016).

The visual representation in Fig. 2 is necessarily simplified, because although the wind structure has azimuthal symmetry the observer also sees imprints on the line profile from the rotational velocity. The relative importance of rotation in determining the line profile shape depends on the location of the line formation region in velocity space (see e.g. Section 4.3). The effect of the rotational velocity on the line profile is discussed explicitly in Section 4.1.1.

4 RESULTS FROM A SIMULATION GRID

We now present results from the full simulation grid, comprising 243 photoionization and radiative transfer simulations with different wind model parameters. The simulations were run on the Cascake Lake nodes of the Cambridge Service for Data Driven Discovery (CSD3) cluster, with an average run-time of ≈200 core-hours per model. We vary the wind mass-loss rate, |$\dot{M}_w$|⁠, the acceleration length R_v, the acceleration exponent, α, the terminal velocity parameter, f_∞ and the inner wind opening angle, θ_min. The outer wind opening angle, θ_max, is adjusted to maintain the desired covering factor Ω/4π (which is set to 0.2 for the initial grid). We chose these specific parameters because they will modify the kinematics and/or line emissivity of the wind by altering the density and velocity along a streamline, with a resulting impact on emission line profiles. These choices are, however, somewhat arbitrary, which should be borne in mind when interpreting our results (see Sections 4.4 and 5.5). For each wind model, we compute synthetic spectra from 5° to 85° at 5° intervals over the wavelength range 500−7500 Å. Convergence in each cell is evaluated by requiring that the electron and radiation temperatures are not changing significantly from cycle to cycle, and the heating and cooling are balanced, both within a relative tolerance of 5  per cent. The convergence criteria are described in the code documentation and by Long & Knigge (2002) and M20. The models generally converged extremely well: the median convergence was 97.6  per cent across the simulation grid. A handful of models did not converge well, particularly for α = 1.5, where the slow acceleration can lead to a dense and optically thick wind base. This structure makes it difficult to achieve good r-packet statistics in shielded regions, with a corresponding minimum convergence of 67.8  per cent. Nevertheless, the median convergence values for each value of α = 0.5, 1, 1.5 were, respectively, 99.2, 97.2, and 96.8  per cent.

4.1 Line profiles and inclination trends

In Fig. 3, we show normalized C iv λ1550 line profiles for a subset (18 of 243) of the simulation grid, colour-coded by inclination. The parameters for each of these labelled models are given in Table A1 in the Appendix. We can see that the spectrum presented in the previous section is not an isolated case. Rather, blueshifted and asymmetric C iv λ1550 line profiles are always observed for some viewing angles in these examples. In tandem with the heuristic demonstration above, this is the first main result of our work: biconical disc winds are successful in creating blueshifted emission lines.

We have focused on the α = 0.5 and α = 1 models in Fig. 3, partly for brevity and partly because the α = 1.5 models are slightly less effective at creating blueshifted asymmetries in the C iv line profile. The α = 0.5 and α = 1 models in Fig. 3 match each other in the other parameters, as can be seen in Table A1. For completeness, the corresponding spectra for α = 1.5 are shown in the Appendix, and the full set of C iv line profiles from all 243 models are given in the Supplementary Material. Inspection of these spectra reveals that blueshifted asymmetries are ubiquitous in spectra that show prominent lines, with the basic conditions for their formation being that the wind has a sufficiently high C iv λ1550 luminosity and that the mid-plane obscures the receding wind cone.

Some interesting trends with inclination can be seen in the line profiles, where we define inclination, θ_i with respect to the symmetry axis, such that a face on disc is seen at θ_i = 0. The lines are narrower at low inclination and broaden towards high inclination. In some cases, at high inclination the blue wing of the line profile can be absorbed, leading to a redshifted asymmetry in the line profile – prominent examples of this behaviour can be seen in panels q and r. As also found by Zamanov et al. (2002) and Gaskell & Goosmann (2016), redshifted line profiles tend to occur at a critical inclination of θ_red ≈ 90° − θ_min, although this is a necessary but not sufficient criterion for forming redshifted asymmetries in our simulations. The reason these redshifted profiles form is actually quite subtle, since simply viewing a receding component to the wind (as occurs when θ_i ≳ θ_red) is not enough to produce a redshifted asymmetry – in fact, without radiative transfer, one would expect a fairly symmetric profile in this case. The redshifted asymmetry is introduced due to the different path-lengths (and therefore different optical depths/absorbing columns) traversed by photons from the blue and red wings of the line. In the receding part of the wind, photons escape to the observer along the normal to their streamline which is roughly the shortest path to the edge of the wind. The blueshifted line photons in the approaching part of the wind have a longer path to travel to escape and so are more readily absorbed, leading to a redshifted asymmetry. Further discussion of the shape of the emission lines and a comparison to the observations can be found in Section 4.2.

To examine the blueshift evolution with viewing angle more concretely, we show C iv blueshift (calculated from equation 1) as a function of inclination in Fig. 4, for the α = 0.5 models (the same plots for the other values of α can be found in Appendix  A). For all models with reasonable strength emission lines, we see a clear evolution from blueshifted profiles at low inclinations and redshifted or symmetric profiles at higher inclinations. Generally speaking, blueshifts are only observed at low inclination, for all wind geometries. To quantify this further, we show the cumulative distribution function (CDF) showing, for each value of α, the fraction of all the model spectra exhibiting blueshifts that are found at an inclination θ < θ_i. We define a ‘blueshifted’ line as having a blueshift exceeding 500 km s⁻¹ for this purpose. A CDF value of 1 means that no blueshifted spectra are found beyond that inclination. For all values of α, at least 70  per cent of the spectra with C iv blueshifts in our model grid are seen at θ_i ≲ 20°, and ≳95  per cent are viewed at θ_i ≲ 40°. Only a handful of blueshifted line profiles are seen at θ_i ≥ 50° for any wind geometry or value of α.

These results demonstrate that blueshifts are preferentially seen at low inclinations, but there are a couple of caveats and complications. First, in these plots we only consider θ_i ≤ θ_min, that is, we have not considered viewing angles looking into the wind cone – since these sightlines often produce BALs (see Section 4.3) – or looking underneath the wind. Along the latter class of sightlines, quite high EW emission lines can be observed as the continuum is often suppressed; we refer the reader to the discussion by Matthews et al. (2016). Our neglect of high inclination viewing angles means that the effect of blueshifts only being produced at low inclination could be artificial. However, this is not the case. Considering only the models with an equatorial disc wind geometry (and so also with a large range of angles with θ_i ≤ θ_min), we still see a clear trend from high C iv blueshift at low inclination to low or negative blueshift at high inclination in the left-hand panel of Fig. 4. In addition, in our adopted geometry it makes intuitive sense that blueshifts would be produce at low inclinations because the mid-plane acts as the obscurer of the red wing. The second complication here is that blueshift can in some cases initially increase with inclination, as occurs in, for example, a number of the θ_min = 45° curves in the left-hand panel of Fig. 4. Again, this is caused by projection effects. As θ_i increases, the outflow velocity component projected into the line of sight can also increase, leading to a slightly stronger blueshift. Neither of these two caveats invalidate our main conclusion here: emission line blueshifts are preferentially produced at low inclination, and this behaviour is an inevitable consequence of the projection effects within, and symmetry axes of, the expected quasar disc wind geometry.

4.1.1 Impact of rotation and line transfer

To assess the impact of rotation and line transfer, we take the two fiducial models, e and n, and, using the same ionization and temperature structure, re-simulate their spectra without rotation or line transfer. Specifically, for the rotation test, we set all rotational velocity components to zero, and for the line transfer test we set the Sobolev optical depth to zero during the actual radiative transfer process (it is still included as a local effect via the Sobolev escape probabilities). The goal of this exercise is not to mimic any physically realistic situation. Instead, the aim is to assess the role rotation and line transfer play in determining the shape of C iv emission lines in our synthetic spectra.

The results from these two tests are shown in Fig. 5. We find that neglecting rotation leads to a single-peaked blueshifted profile which is also narrower. The impact is more pronounced for Model e because the line is formed closer to the wind base where the kinematics are more dominated by rotation. In contrast, rotational velocity has a fairly modest effect on the model n line profile, because outflow is more important in the line formation region. The overall line profile can be thought of as a convolution of a rotational component with an outflowing one, with the relative importance of each determined by the kinematics at the point where C iv λ1550 is emitting efficiently – although, in reality, the line formation region is stratified in height and thus velocity, as discussed in Section 4.3.

The neglect of line opacity within the flow leads to a modest, velocity-dependent suppression of the overall line flux. This suggests there is a resonant scattering component to the line profile, caused by scattering of C iv line photons into the line of sight. In model e, this occurs preferentially at blueshifted velocities and is responsible for a horn-like feature, reminiscent of one side of a double-peaked disc emission line profile. In model n, the scattering contribution is less confined to one part of the line profile although it does occur exclusively at velocities ≲1000 km s⁻¹. These results suggest that radiative transfer effects within the wind are important for determining the overall line profile, and so only knowing the emissivity profile, say, of emission lines within the wind is not sufficient for an accurate characterization of their shapes.

4.2 Comparison to observations

4.2.1 Skew and line profile shape

One characteristic feature of the line profiles presented in Fig. 3 is that the peak in the line profile tends to be skewed towards the blue, such that the velocity of the peak is more highly blueshifted than the ‘blueshift’ calculated from equation (1). Using skew, as defined by equation (2), we can quantify this effect for both the observations and models, the results of which are presented in Fig. 6 for the α = 0.5 and α = 1 models. In the composite spectra from R20, we find a roughly linear relationship between skew and blueshift, with more blueshifted lines having more negative skews. In other words, and as can be seen from Fig. 1, weak line quasars with strong blueshifts tend to have strong asymmetries and a more extended blue wing, with their line peaks lying redwards of the wavelength that bisects the line flux. This is generally the opposite to what we see in the models, where high-blueshift spectra tend to have positive skews in Fig. 6. Thus, while we are successful in producing blueshifts of a couple of thousand km s⁻¹ in our models, the line profiles are normally skewed in the wrong direction compared to the observations. Furthermore, in the models, there is a clear weak positive correlation between blueshift and skew, in the opposite sense to the anti-correlation in the composite spectra. This result might hint at problems with the wind geometry or kinematics in our model compared to real outflows in quasars. However, it is also important to stress once again that Fig. 6 is in some sense an ‘apples versus oranges’ comparison, since we have simulated a quasar with fixed accretion rate and mass, and we are comparing to composite spectra built from a population of quasars spanning a range of masses and Eddington fractions.

There are a few other effects in Fig. 6 that are worth discussing. The first can already be inferred from Fig. 4, that the most extreme blueshifts in the model do not exceed 2000 km s⁻¹, whereas the composite spectra (and the overall C iv emission space) extends to ∼3000 km s⁻¹. Our wind model cannot reproduce these extreme blueshifts, suggesting that these outflows are particularly dense or powerful such that C iv line formation can take place efficiently at high velocity (see discussion in subsequent subsection). This apparent discrepancy is perhaps expected. The ∼3000 km s⁻¹ C iv blueshifts are observed only in a small population of quasars with Eddington fractions of ≳0.5 and slightly higher masses than we have considered here (Temple et al. 2023). These quasars would be expected to have softer SEDs (or higher UV to X-ray flux ratios) more conducive to line-driving, which, combined with the high Eddington ratios, might be expected to lead to more powerful and/or faster outflows (e.g. Giustini & Proga 2019). The second notable point is that there are points more or less overlying some of the observational data. However, inspection of these models reveal that none of them produce line profiles with similar shapes to the observations, as can be seen from the full set of line profiles in the Supplementary Material. Finally, we note that different wind geometries – denoted by different colours in Fig. 6 – produce different patterns in this parameter space, with more polar wind geometries tending to have more positive skews. The fact that there is a relationship between the line profile shape and wind geometry merits further investigation, implying that different or more complex wind geometries to those considered here could lead to a better agreement with observed spectra (see Section 5.1).

4.2.2 Equivalent width and line luminosity

When inspecting all 243 of the C iv λ1550 line profiles from our simulation grid (see Supplementary Material), we find that blueshifts are always produced at low inclination as long as the emission line is obviously present. In a number of cases, the wind is over- or underionized, and the continuum is more or less featureless in the UV. Thus, while Fig. 3 shows a favourable selection of models that all produce strong lines, the blueshift formation is an inevitable consequence of obscuring the receding part of the wind – as long as the lines are strong enough and the obscuring mid-plane is present, a blueshifted asymmetry is formed. However, ensuring enough of the wind lies in the right ionization state to form emission lines of sufficiently high EWs (that is, comparable to those in observed spectra) can be a problem, as we have discussed in previous work (Matthews et al. 2016, M20).

As we have described and illustrated in Section 2.1, there is a seemingly profound relationship between EW and blueshift in quasar spectra. The true C iv emission space is built up from quasars with a range of masses and bolometric luminosities, whereas our models are run for one specific value of each. As a result, it would be a little misleading to compare directly the range of EWs and blueshifts in our models to those from a population of quasars, since our intention is not to simulate such a population. Instead, we can compare the overall luminosities and EWs of C iv λ1550 computed from our synthetic spectra to observations. This exercise has already been carried out with similar models in our previous work (M20, see their fig. 6) where we found that models that avoided overionization produce line luminosities and EWs similar to those in observed spectra of quasars of comparable luminosity. The same is true with our simulation grid here: there are still a substantial number of models that do not produce noticeable emission lines, because they are either under- or overionized. Considering all models, the median C iv line luminosity at viewing angles ≤θ_min is L₁₅₅₀ ≈ 1.4 × 10⁴⁴ erg s⁻¹, spanning a wide range of luminosities (≈10⁴²−10^45.5 erg s⁻¹). For the subset of models shown in Fig. 3, the line luminosities at the same viewing angles are a little higher, since these models have been selected to have obvious C iv line profiles, with a median of L₁₅₅₀ ≈ 5.4 × 10⁴⁴ erg s⁻¹. For reference, SDSS quasars with λL₁₃₅₀ ≈ 10⁴⁶ erg s⁻¹ typically have C iv line luminosities L₁₅₅₀ ≈ 3 × 10⁴⁴ erg s⁻¹, similar to the median line luminosity in our models.

4.3 Velocity structure and BAL formation

To allow an inspection of the kinematics of the line formation region, we show a ‘phase space’ plot of vertical velocity, v_z, and electron density, n_e, in Fig. 7, with the panels matching those in Fig. 3. Each point shows the location in n_e, v_z parameter space for each cell in the wind model, colour-coded by the local C iv λ1550 line luminosity in each cell. This quantity includes the local Sobolev escape probability but the line photons can be absorbed further out in the wind, so do not necessarily escape to infinity. The bright regions in the plot show the typical velocities and densities where the C iv λ1550 line is forming. There is some variation in the velocity and density structure across the grid (due to changes in the grid parameters θ_min, |$\dot{M}_w$|⁠, R_v and f_∞), but typically the main line formation region spans the density range 8.5 ≲ log (n_e) ≲ 10.5 and the velocity range 300−3000 km s⁻¹. To make things more quantitative, we can define v₉₀ as the v_z such that 90  per cent of the C iv line emission comes from regions with v_z < v₉₀. In all models shown in Fig. 7, v₉₀ < 3050 km s⁻¹, and the mean and median v₉₀ values are 1311 and 1335 km s⁻¹, respectively.

The location of the line formation region in the n_e, v_z plane makes clear why characteristic blueshift velocities of a few 1000 km s⁻¹ are observed in the spectra computed from these wind models; at higher velocities the density is lower and the emission-line emissivity is much lower. The line emissivity is driven by two factors. The first is ionization state – at low densities and high velocities the wind is overionized (mostly in C v and upwards) and C iv line formation is ineffective. The second is density or emission measure. Since the line is collisionally excited, its emissivity is proportional to |$n_e^2$|⁠, so even for fixed ionization state a drop in density leads to a corresponding drop in emissivity. Both these factors act to confine the C iv line formation region to near the base of the outflow where velocities are relatively modest and densities are relatively high. The line formation region can be moved around in the flow by, for example, modifying the velocity law (see Section 4.4).

The wind does eventually accelerate towards its terminal velocity, becoming a more tenuous low-density plasma with velocities exceeding 10 000 km s⁻¹. We find, for observers looking directly into the wind cone, that C iv BALs are often formed in the extended high-velocity flow. In Figs 8 and 9, we show the simultaneous formation of C iv BALs and blueshifts in the same outflow for Models e and n. In both cases, spectra in the rest-frame UV from low to high inclination are shown in the right-hand panel with offsets applied for clarity, and the viewing angles through, and ionization structure of, the wind are shown in the left-hand panel. At low inclination, blueshifted emission lines are observed, and at angles θ_min ≲ θ_i ≲ θ_max BALs are observed in C iv λ1550 and N v λ1240, and often also S iv λ1400. As we found previously (Matthews et al. 2016), LoBAL features in, e.g. Al iii λλ1855,1863 and even Fe ii and Fe iii absorption lines can be observed at high inclination angles that look through the low-ionization, shielded part of the wind.

Consistent with our earlier work (Higginbottom et al. 2013; Matthews et al. 2016), BALs are only formed when the ionization state is sufficiently moderate and a significant line-opacity is present in the flow. As such, BALs typically form in the same models that are successful in producing strong emission lines, although the BAL formation region is significantly more extended than the emission line region. Figs 8 and 9 represent another key result of our work: BALs and C iv-emission blueshifts can form in the same wind structure within a geometric unification scheme. This result is discussed critically in the context of observational results in Section 5.2.

4.4 Dependence on model parameters: the velocity law is key

Our simulation grid is built from five main model parameters, and the synthetic spectrum is determined by the combination of these plus observer inclination, θ_i. We have already discussed the impact of θ_min and θ_i. To investigate the sensitivity of the C iv line profile shape to the model parameters, we take Model e as a starting point and then take steps away from this model in the four remaining parameters (α, R_v, f_∞, |$\dot{M}_w$|⁠). In addition, we also varied the obscuring disc size, R_mid and the covering factor, Ω/4π. To provide slightly finer resolution, we supplement our main simulation grid with a few additional model runs so as to have spectra at five values of each parameter.

The results of this ‘step’ parameter search are shown in Fig. 10, designed to give the reader a feel for how the C iv emission line varies with each of these quantities. From this figure, we can see varying degrees of impact on the line profile shape and line EW. For example, decreasing the acceleration length (or increasing the terminal velocity through the parameter f_∞) acts to weaken the emission line but has little impact on the line profile shape. By contrast, the acceleration exponent has a significant effect, changing the peak velocity of the line as well as its skew and EW. This dependence of the line profile shape on the acceleration exponent, and therefore, the shape of the velocity law is particularly interesting, so we discuss it further in Section 5.1. Additionally, we can see that the extent of the opaque mid-plane, R_mid, is important in determining whether blueshifts form, as expected from the discussion in Section 3. In particular, blueshifts start to emerge once R_mid ≳ 10¹⁸ cm, a threshold value matching the characteristic horizontal distance of C iv λ1550 line formation in this model; we discuss the opaque mid-plane and this physical scale further in Section 5.3.

The reason for the change in line profile shape with the acceleration exponent, α, can be understood more clearly by examining how changing α affects the location of the line formation region in velocity space. In Fig. 11, we show the C iv line profiles for each of the five values of α (top panels), together with the location of the wind cells in v_l, v_ϕ parameter space, colour-coded by line luminosity as in Fig. 7 (bottom panels). Additionally, we show the poloidal velocity along a streamline from equation (3) for each value of α in the right-hand panel, as well as the rotational velocity along the streamline as calculated from conservation of angular momentum. As α is increased, the wind accelerates more slowly and the line formation region moves to higher ratios of v_l/v_ϕ, producing emission-line profiles that are narrower, single-peaked, and more dominated by outflow. By contrast, for α = 0.5 the line formation is mostly confined to regions close to the disc where v_ϕ is high.

Further intuition can be obtained here by considering a simple toy model where the C iv line forms at a fixed ionization parameter ξ. Assuming that changing the velocity law does not alter the line formation radius all that much, and considering fixed ionizing luminosity and wind mass-loss rate, then the ionization parameter is set by the density. The density, in turn, is largely determined by the poloidal velocity of the wind (Shlosman & Vitello 1993). In this case, one might expect that the line formation region occurs around a given fixed value of v_l. The shape of the line profile, and the relative dominance of rotation and outflow components, is then set by v_l/v_ϕ. Increasing α means that the wind accelerates more slowly, but the variation of v_ϕ along a streamline is fixed (see right-hand panel of Fig. 11). Thus, increasing α causes the ‘critical’ value of v_l to be reached at lower values of v_l/v_ϕ, and thus the line profile becomes more single-peaked and less influenced by rotational effects.

4.5 Impact of SED shape and anisotropy

For our initial simulation grid and fiducial model, we used the same isotropic ‘disc + power-law’ SED as M20. The anisotropy and shape of the illuminating SED is likely to play an important role in determining the observed spectrum: the flux anisotropy changes wind dynamics (Dyda & Proga 2018; Williamson, Hönig & Venanzi 2019) and emission line EWs (Risaliti, Salvati & Marconi 2011; Matthews et al. 2017), while the frequency distribution modifies both the force multiplier (Dannen et al. 2019) and the observed emission line ratios (e.g. Leighly 2004; Richards et al. 2011).

In this work, we cannot account for dynamic effects directly, but we can investigate the sensitivity of our results to SED changes for a fixed wind structure. To do this, we take models e and n, and fully re-simulate the ionization structure and resulting C iv λ1550 emission-line profile for four additional SED choices: an identical SED shape with an anisotropic accretion disc (following Matthews et al. 2016), an isotropic SED based on the Richards et al. (2006) template, and two SEDs computed with the qsosed model of Kubota & Done (2018). The qsosed models were taken from the Temple et al. (2023) simulation grid and are described by mass, Eddington-normalized accretion rate (⁠|$\dot{m}$|⁠) and dimensionless spin (a_*). We pick two models: qsosed1, with a BH mass of 10^9.5 M_⊙, |$\dot{m}=0.468$| and a_* = 0, and qsosed2, with a BH mass of 10⁹ M_⊙, |$\dot{m}=0.15$| and a_* = 0.998. These two model spectra are quite different in terms of the hardness or colour of the SED and produce very different fluxes. We renormalized all of our model SEDs to have identical bolometric luminosities to that of our original model. As a result, the SEDs used are not necessarily physically consistent or realistic, but rather are designed to illustrate the variation of the line profiles for a wide range of spectral shapes and ionizing fluxes.

The results from our SED investigation are shown in Fig. 12. In the right-hand panel, we show the four SEDs, while in the left-hand panels we show how the line profiles change in each case, for both models e and n. Most of the SED changes have a relatively modest impact on the C iv line profile, with the results from the Richards et al. (2006) and qsosed1 SEDs showing rather similar profile shapes and line EWs to the fiducial model case. The main effect of the anisotropic disc is to increase the continuum flux (due to foreshortening and limb darkening), which acts to decrease the line EW slightly. Easily the most dramatic change is observed when using the rather extreme qsosed2 SED, which produces a significantly stronger line in the model n case. This SED has a significantly stronger EUV component as well as a much weaker continuum in the near UV. The key frequencies corresponding to the line wavelength and the C iii→C iv and C iv→C v ionization edges, as well as the C iv inner shell edge, are marked in the right-hand panel of Fig. 12. The weaker continuum underlying the line means the EW is significantly higher, but there is also an effect on the line luminosity, which increases. In tandem, the shape of the line profile has changed and now has a slight redshifted asymmetry due to absorption of the blue wing by the wind. However, these more dramatic changes are only seen for this drastically different SED shape, which has much higher fluxes around the key ionization edges (see right-hand panel of Fig. 12). This SED is unlikely to be particularly realistic for a typical quasar,

Overall, we can draw two main conclusions from our SED investigation. First, that significant changes in SED shape, as might be seen across a population of quasars, can clearly alter the C iv blueshift and EW; this result is interesting given the apparent dependence of on C iv emission space location on SED shape (Richards et al. 2011; Rankine et al. 2020; Temple et al. 2023). Secondly, that our overall results are not particularly sensitive to our adopted SED, since three reasonable choices with the same bolometric luminosity but different soft X-ray excesses and ionizing luminosities all produce blueshifted C iv emission lines with similar shapes and EWs.

5 DISCUSSION

5.1 Line profiles as a probe of disc wind physics

We found in Section 4.4 that the velocity law – in particular the value of α – can have a dramatic impact on the observed shape of the C iv λ1550 line profile, and we discussed how the line formation region moves around in velocity space with these changes in α. In a general sense, this finding suggests that emission-line profile shapes could be used as a probe of the velocity law of the wind. Although the detailed wind velocity law from different wind driving mechanisms is not well known, there are generic reasons to expect them to have quite varied kinematic structure, as different driving forces should produce winds with correspondingly different acceleration scale lengths and terminal velocities. Additionally, magnetocentrifugal winds conserve angular velocity along streamlines, as opposed to specific angular momentum as in the case of thermal or radiation driven outflows. Thus, magnetocentrifugal winds would be expected to have higher values of v_ϕ/v_l in the line formation region compared to other driving mechanisms (all other parameters being equal). Further modelling is needed here, but this argument clearly demonstrates the link between line profile shape and wind driving mechanism. In fact, the presence of outflow-dominated kinematic components in quasar spectra – necessitating a line formation region where v_l > v_ϕ – could already be pointing towards, say, line-driven winds as the origin of these emission lines.

One specific aspect of the line profile shape that changes with varying α (and also R_v and |$\dot{M}_w$|⁠; see Fig. 10) is whether or not the line is single peaked. Our work here demonstrates that disc winds can produce both double- and single-peaked lines, with the observed profile depending on both inclination and the wind parameters. The basic mechanism for the single-peaked behaviour is, once again, the formation of the emission line in an outflow dominated region where, approximately, v_l > v_ϕ, although the details also depend on the projected velocity and resonant line scattering. We found similar behaviour in our modelling of disc winds in accreting white dwarfs. There the line formation region moved up in the flow as R_v was increased and could cause a transition from double- to single-peaked Balmer lines. This behaviour is distinct from the oft-discussed mechanism proposed by Murray & Chiang (1996, 1997), in which radial velocity gradients modify the escape probabilities and lead to a single-peaked profile.

The fact that line profile shape depends so clearly on the wind kinematics rather begs the question: what kind of outflows could better mimic the observed broad emission-line shapes? We cannot conclusively answer this question at this stage, but it would seem reasonable that alternative prescriptions for quasar outflows could produce a peak of emission that is redwards of the blueshift velocity. It would be interesting in future work to consider outflows with more complex kinematics and multiphase structures. One example is the ‘quasar rain’ model suggested by Elvis (2017). Alternatively, there are various ‘failed wind’ scenarios that have been explored in both dust-driven (Czerny & Hryniewicz 2011; Galianni & Horne 2013; Czerny et al. 2017; Naddaf & Czerny 2022; Naddaf et al. 2023) and line-driven (Proga & Kallman 2004; Giustini & Proga 2019) flows, which are also relevant to the two-component BLR discussed in the next subsection. These more complex or clumpy wind geometries can represent challenges to the MCRT methods we have used, but are clearly worthy of further investigation.

5.2 C iv emission space, the BAL connection, and the BLR

In Section 4.3, we showed that BALs and blueshifted emission lines could be produced by the same disc wind structure, which is interesting in the context of the results of R20. Using an emission-line reconstruction scheme, R20 found that BALs and non-BALs overlap across much of C iv emission space and also demonstrated that the BAL troughs with the highest maximum velocities tended be found at high C iv blueshifts. The idea that BALs and emission lines are formed in the same disc wind structure has been explored extensively in various ‘unification’ schemes and modelling efforts (e.g. Murray et al. 1995; Elvis 2000; Matthews et al. 2016), but our work explicitly demonstrates that MCRT simulations can simultaneously produce not just emission lines and BALs, but also the blueshifted asymmetry of the emission lines. However, there are various caveats and complexities to consider in addition to the skew problem already discussed in Section 4.2.1. Inspection of Figs 8 and 9 reveals that the line profiles close to BAL viewing angles look quite different to those at low inclinations. Thus, just because BALs and blueshifts are formed in the same outflow does not mean those BAL and non-BAL systems will have similar observed emission line spectra, because of the inclination dependence of, e.g. line widths, EWs and blueshifts. A possible solution here is to consider a more filamentary or spoke-like structure for the wind model, as discussed by Yong et al. (2018).

In our work, we have focused almost exclusively on the C iv λ1550 emission line. This is partly for simplicity, and partly because the blueshifted asymmetry in this line is the cleanest, strongest and best-studied signature of its type in quasar spectra. We presented full UV spectra for two models in Figs 8 and 9, in which various UV lines such as Ly α, N v λ1240, He ii λ1640, S iv λ1400, and C iii] λ1909 can be clearly seen, but we have not discussed the kinematic signatures in these lines in detail. Such signatures are astrophysically relevant – for example, Temple et al. (2021) show that, in quasars with high C iv blueshifts, the Si iv, Ly α, and N v lines all have blueshifted asymmetric structures resembling C iv. In addition, the C iv blueshift is often defined with respect to a line like Mg ii λ2800  that typically has a centroid close to the systemic velocity. A detailed investigation of the range of emission lines profile shapes in wind models would be an interesting avenue for future work. Nevertheless, it does seem fairly clear from the line ratios and profile shapes in Figs 8 and 9 that it is quite challenging to get a biconical wind mode of the type considered here to reproduce the full spectrum of a high-blueshift ‘windy’ quasar, even if disc winds can imitate ‘locally optimally emitting’ behaviour (Baldwin et al. 1995) by spanning a range of densities and ionization states, as suggested by M20.

It is probably more reasonable to think of the wind-formed emission lines simulated here as corresponding to one component of a multifaceted broad-line region. If we consider the simplest possible multicomponent BLR, the degree of blueshifted asymmetry (and the location in C iv emission space) would be partly driven by the relative strength of the ‘windy’ and ‘peaky’ kinematic components (see e.g. Temple et al. 2021, 2023). The physical origin of these two components could then correspond to a line- or magnetically driven disc wind, as considered here, and a classic BLR-like component – perhaps a rotationally dominated turbulent plasma, failed wind or population of clouds – that produces a symmetric line. This two-component BLR is clearly a toy model; it seems more likely that the BLR is a stratified structure, with complex, probably turbulent, kinematics. Further modelling combined with novel data analysis (for example, attempts to decompose the kinematic components) will thus be essential to reveal its true nature.

5.3 The opaque mid-plane

One of the key ingredients for forming blueshifts in our modelling is the presence of an opaque mid-plane out to large radii, such that R_mid exceeds the characteristic radius of line formation. It was noted by Wilkes (1984) nearly four decades ago that the obscuration of the red wing could be due to ‘...the presence of a large disc of gas and dust, probably an extension of gas accreting on to a central black hole,’ extending out to ∼10¹⁸ cm (∼1/3 pc). Remarkably, this is the same value of R_mid at which blueshifts start to form in the top left-hand panel of Fig. 10, as dictated by the location of emission line formation. The critical value of R_mid likely depends somewhat on the specific model, since we know from Fig. 11 that the line formation region moves around in the wind (in both physical and velocity space) for different wind parameters. Nevertheless, it is interesting to compare this radius to some characteristic radii in the quasar system: for example, taking a 10⁹ M_⊙ black hole accreting at |$\sim 20~{{\ \rm per\ cent}}$| of Eddington, the dust sublimation radius is R_dust ∼ 1 pc (Barvainis 1987) and the self-gravitation radius is R_sg ∼ 0.02 pc (King 2016). Thus, R_sg is significantly less than, and R_dust is broadly comparable to, the required opaque mid-plane extent. This result is supported observationally by the size of the BLR being similar to R_dust in many AGN including high mass quasars 3C 273 and IRAS 13349 + 2438 (Gandhi, Hönig & Kishimoto 2015). We note, however, that it is not clear whether the BLR radius estimated from reverberation mapping is relevant to the line formation radius for the C iv blueshifted component in quasars. Indeed, Temple et al. (2023) argue, based on a need for high ionizing fluxes, that the ‘windy’ component of a two-component quasar BLR is formed much closer to the central black hole.

If the line formation radius of the blueshifted C iv is comparable to that in this work, it is fairly likely that the obscurer is a self-gravitating, dusty disc or vertically extended torus, rather than a geometrically thin, optically thick accretion disc. In this case, the obscuring mid-plane may also be playing a role in forming the blueshifts sometimes observed in the narrow lines. In addition, a dusty midplane would have a frequency-dependent opacity. Alternatively, it may be that dust opacity within or surrounding the wind itself has a role to play in absorbing the red wing of broad emission lines, particularly if the wind is dust-opacity driven. Thus, while the presence of an extended obscuring disc is a fairly natural consequence of the axisymmetric quasar system, future work might focus on incorporating dust radiative transfer through a physically motivated model for a dusty disc (and/or dusty outflow), rather than simply treating the mid-plane as a ‘brick wall’.

5.4 Interpretation of line asymmetries

We close this section by making some final, rather general comments on the interpretation of asymmetries in emission lines, which are relevant to a host of astrophysical settings. Blueshifted asymmetries in emission lines are not smoking guns for outflows – they tell us that the line-emitting material is moving towards us, but, unlike BALs, they do not tell us whether that material is located between the observer and the continuum source. Indeed, in our work the disc wind was able to produce both blue- and redshifted emission line signatures from the same wind geometry, highlighting the potential for confusion.

Velocity-resolved reverberation mapping results provide important constraints on the kinematics of the BLR (e.g. Peterson 1993; Denney et al. 2009; Pancoast, Brewer & Treu 2011; Waters et al. 2016; Mangham et al. 2017; De Rosa et al. 2018), although most studies focus on local, lower luminosity AGN, which are very different systems to the high-blueshift quasar population. Nevertheless, results can be informative for our work; for example, in a recent study, Oknyansky et al. (2021) find spectacular asymmetric lines during a changing-look event in the Seyfert galaxy NGC 3516. At first glance, the asymmetries in H α and H β they observe have very similar shapes to some of the blueshifts we have found in our models. It is therefore tempting to interpret those asymmetries as coming from a disc wind driven during the changing-look event. However, the reverberation lags in the lines have a red-leads-blue signature which is characteristic of inflow (Welsh & Horne 1991). The picture is complicated even further, because recombination lines in particular can have complex responses and transfer functions (Korista & Goad 2004; Mangham et al. 2017, 2019) which can even lead to red-leads-blue signatures from outflows, as shown by Mangham et al. (2017, see also Waters et al. 2016; Yong et al. 2017). The Oknyansky et al. (2021) results are in some sense cautionary, since the NGC 3516 data show one signature (emission line blueshift) commonly interpreted as evidence for an outflow, and one signature (red-leads-blue reverberation) commonly interpreted as evidence for an inflow.

Given all the above, one might legitimately ask: are blueshifted emission lines telling us anything at all about accretion disc winds? While it is true they cannot be considered ‘smoking guns’ of outflowing material in the same way as blueshifted absorption lines, there are nevertheless at least three pieces of evidence that point towards their formation in outflows. The first is simply physical plausibility; we and others (Richards et al. 2002; Chajet & Hall 2013, 2017; Yong et al. 2017; GRAVITY Collaboration et al. 2020) have demonstrated that a reasonable disc wind geometry can produce blueshifted asymmetries. We do not know of an alternative, plausible, physically motivated model for their formation² (e.g. in a plasma whose kinematics are dominated by rotation or inflow). The second piece of evidence is the behaviour of BAL properties in C iv emission space. R20 find that the most extreme BALs (i.e. highest Balnicity index [BI] and largest maximum trough velocities) are typically associated with the largest C iv blueshifts when the spectrum is reconstructed using mean-field independent component analysis. While this relationship does not necessitate a direct physical association, the simplest explanation is that C iv blueshifts are created by the BAL outflows themselves. Finally, there is the general point that, from the prevalence of C iv BALs, we know quasars produce outflows carrying a significant C iv ion population, so it seems rather natural to associate these outflows with the C iv emission lines. This argument is strengthened by our work, since we have shown that C iv BALs and emission line blueshifts can be formed in the same disc wind under reasonable assumptions (although see also discussion in Section 5.2).

5.5 Limitations and model assumptions

Our work has a number of limitations, which we briefly discuss as well as referring the reader to our earlier work (Matthews et al. 2016, M20).

Given the obvious impact of wind geometry and the velocity law on the line profile shape, it seems clear that adopting the Shlosman & Vitello (1993) wind parametrization is a limitation of this work. In particular, it may be possible to get a better match to observed quasar spectra with a different wind model. It would be valuable to test model prescriptions that are more physically motivated, in the sense that the wind streamlines and opening angles are obtained from a self-consistent solution of the wind dynamics. Similarly, although our choices of variable wind parameters and the values of the grid points in our simulation grid are sensible, there are also feasibly other model parameters or regions of parameter space that should be investigated.

In terms of the physics in the simulations, the main limitations are probably the treatment of clumping and the SED. Our ‘microclumping’ treatment, in which inhomogeneities in the flow are implemented through a single volume filling factor, is necessarily simplistic and clearly does not capture the physics of multiphase, clumpy outflows. In terms of the SED, we have assumed an isotropic SED, and while we have investigated the sensivitity of the C iv line profile to both the SED shape and anisotropy (Section 4.5), the treatment of the continuum radiation field can dramatically affect the ionization state of the wind and the EW of the lines in the synthetic spectra.

6 CONCLUSIONS

We have demonstrated how blueshifts can form in biconical disc winds (e.g. Figs 2 and 3), and used MCRT and photoionization simulations to examine trends with inclination, wind geometry, and other wind parameters. Beyond this result, our main conclusions are as follows:

• We find that a key condition for blueshift formation is that the disc mid-plane is opaque out to radii beyond the line formation region, so that the receding wind bicone is obscured (Section 3; Figs 2, 3, and 10).

• Blueshifts are only produced at relatively low (that is, more face on) inclinations, with ≈70  per cent of the spectra with C iv blueshifts seen at θ_i ≲ 20°, and ≈95  per cent at θ_i ≲ 40° (Section 4.1; Fig. 4).

• We find that a combination of ionization and emission measure effects act to confine the C iv line formation region to near the base of the outflow (Section 4.3; Fig. 7), where velocities are typically 300−3000 km s⁻¹. Further out in the flow, densities are too low, and the C iv λ1550 emissivity drops dramatically. This mechanism provides a natural explanation for the size of the observed C iv blueshifts and shows that emission line velocities do not necessarily provide any information about the terminal velocity in the flow. It also suggests that higher C iv blueshifts would require higher mass-loss rates or lower volume filling factors, such that the lines could form at higher velocity.

• We demonstrate that BALs and emission-line blueshifts can be formed in the same outflow (Section 4.3; Figs 8 and 9), with BALs observed at higher inclinations when the viewing angle intersects a significant portion of the wind cone (θ_i ≳ θ_min), and emission line blueshifts at low inclinations (θ_i ≲ θ_min). This result has interesting implications for our understanding of the BAL versus non-BAL dichotomy and, in particular, the C iv emission space. However, it is challenging to construct a workable biconical model in which the emission-line properties in the BAL and non-BAL quasars are as similar as in the observational data (R20).

• Despite our success in reproducing blueshifted emission lines within a disc wind model, we cannot match the detailed line profile shapes (Section 4.2.1). In particular, the model spectra tend to be ‘positively skewed’ in the sense that the line peak lies bluewards of the blueshift calculated from the wavelength bisecting the line flux; observed spectra display the opposite behaviour. This discrepancy may suggest a fundamental problem with a disc wind origin for C iv blueshifts, but it is probably more likely that our biconical model is limited in its parametrization and a poor imitation of a real, dynamic quasar wind.

• Redshifted emission lines can also be produced when the wind region on the far side of the black hole is viewed at such an angle that the plasma is receding from the observer (Section 4.1). This starts to occur above a critical viewing angle of θ_red ≈ 90° − θ_min, a necessary but not sufficient criterion for redshifted line emission. Radiative transfer effects (i.e. preferential absorption of the blue wind) within the wind are critical for forming these asymmetries even when θ_i ≳ θ_red.

• We find that the shape of the velocity law has a significant impact on the shape of the C iv line profile, which we demonstrate through variation of the acceleration exponent α (Section 4.4; Figs 10 and 11). The velocity law has a twofold effect on the line profile by changing the density structure, and therefore location of the line forming region in velocity space, and also altering the projected velocity into the line of sight. A critical factor is how quickly v_l increases along streamlines with respect to the decrease of v_ϕ. Importantly, the impact of the wind velocity profile suggests that emission-line blueshifts and asymmetries can be used as a probe of disc wind driving mechanisms in quasars and other accreting systems.

Future work should move beyond the limited parametrization of a disc wind we adopted, and, in an ideal world, construct ab initio models with different driving physics. When confronted with observations, such an approach will help establish how blueshifted emission lines are formed and more generally help understand the physical mechanisms driving the diversity of line profile shapes and associated BLR kinematics. Regardless of the driving mechanism, accurate radiative transfer calculations will be critical in order to produce consistent synthetic spectra. In addition, our work suggests that the emission lines can in principle be a powerful probe of the underlying physics of outflows/winds in quasars and other accreting systems.

ACKNOWLEDGEMENTS

We thank the referee for reviewing our paper and their constructive comments. We would like to thank Andrew Sellek, Stuart Sim, Nico Scepi, and Teo Muñoz-Darias for helpful discussions. JM acknowledges support from a Royal Society University Research Fellowship and, previously, the Herchel Smith fund. ALR acknowledges support from a UKRI Future Leaders Fellowship (grant MR/T020989/1). MJT acknowledges support from a FONDECYT postdoctoral fellowship (3220516). This work was performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (www.csd3.cam.ac.uk), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP/T022159/1), and DiRAC funding from the Science and Technology Facilities Council (www.dirac.ac.uk). We gratefully acknowledge the use of the following software packages: OpenMPI (Gabriel et al. 2004), matplotlib (Hunter 2007), and astropy (Astropy Collaboration 2013, 2018). For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

Funding for the SDSS IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss4.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics| Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

DATA AVAILABILITY

The bulk of the simulation data used in this article is public in a github repository at https://github.com/jhmatthews/blueshifts, with an associated Zenodo record and digital object identifier (Matthews 2023). The observational data underlying this article come from the SDSS Data Release 16. The composite quasar spectra are available in the Supporting Information included with Stepney et al. (2023). Any other derived data products are available from the authors upon reasonable request. The Python radiative transfer code is an open-source project available at https://github.com/agnwinds/python.

Footnotes

References

APPENDIX: MODEL PARAMETERS AND ADDITIONAL PLOTS

In Table A1, we present the model parameters for the 18 members of the simulation grid shown in Figs 3 and 7, together with labels matching those in the figure panels. In addition, we show the nine models with α = 1.5 with the other parameters matching the corresponding α = 0.5 and α = 1 models. Fig. A1 shows the line profiles (left-hand panel) and wind structure in n_H − v_z space (right-hand panel) from these nine models with α = 1.5. Finally, in Fig. A2, we show C iv blueshift as a function of inclination for the α = 1 and α = 1.5 models, analogously to the left-hand panel of Fig. 4.