{'fraction_non_alphanumeric': 0.031690539775026304, 'fraction_numerical': 0.022351703487901595, 'mean_word_length': 4.532909465389093, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Outflow and hot dust emission in broad absorption line quasars

We have investigated a sample of 2099 broad absorption line (BAL) quasars with z=1.7-2.2 built from the Sloan Digital Sky Survey Data Release Seven and the Wide-field Infrared Survey. This sample is collected from two BAL quasar samples in the literature, and refined by our new algorithm. Correlations of outflow velocity and strength with hot dust indicator (beta_NIR) and other quasar physical parameters, such as Eddington ratio, luminosity and UV continuum slope, are explored in order to figure out which parameters drive outflows. Here beta_NIR is the near-infrared continuum slope, a good indicator of the amount of hot dust emission relative to accretion disk emission. We confirm previous findings that outflow properties moderately or weakly depends on Eddington ratio, UV slope and luminosity. For the first time, we report moderate and significant correlations of outflow strength and velocity with beta_NIR in BAL quasars. It is consistent with the behavior of blueshifted broad emission lines in non-BAL quasars. The statistical analysis and composite spectra study both reveal that outflow strength and velocity are more strongly correlated with beta_NIR than Eddington ratio, luminosity and UV slope. In particular, the composites show that the entire C IV absorption profile shifts blueward and broadens as beta_NIR increases, while Eddington ratio and UV slope only affect the high and low velocity part of outflows, respectively. We discuss several potential processes and suggest that dusty outflow scenario, i.e. dust is intrinsic to outflows and may contribute to the outflow acceleration, is most likely. The BAL quasar catalog is available from the authors upon request.


INTRODUCTION
Outflows appear to be a common phenomenon in quasars. Strong outflows may carry away huge amounts of material, energy and angular momentum and are believed to be one of the most important feedback processes connecting their central engines and host galaxies (e.g. Hopkins & Hernquist 2006;Hopinks & Elvis 2010). For instance, outflows are thought to be able to regulate the growth of supermassive black holes (SMBH) and star formation in host galaxies (see Antonuccio-Delogu & Silk 2010 for a recent review), chemically enrich the interstellar medium (ISM) of host galaxies and the surrounding intergalactic medium (IGM) (e.g., Collin & Zahn 1999;Veilleux et al. 2005). Lots of studies tried to closely tie up outflows with the quasar fundamental parameters, such as Eddington ratio (Boroson 2002;Gangly et al. 2007;Zhang et al. 2010;Wang et al. 2011;Marziani & Sulentic 2012), the intrinsic spectral energy distribution (SED, e.g. Leighly & Moore 2004;Fan et al. 2009;Richards et al. 2011;Baskin et al. 2013) and even gas metallicity (e.g. Wang et al. 2012). Therefore outflows are important for our understanding of the evolution of the central SMBHs and the connection with host galaxies.
Quasar outflows manifest themselves in blueshifted ab-sorption and emission, such as broad absorption lines (BALs; Weymann et al. 1991) and blueshifted broad emission lines (BELs, Gaskell 1982). It is well known that BAL quasars in general have redder ultraviolet (UV) continua than non-BAL quasars, and quasars with BALs from low-ionization state species (e.g. Mg II and Fe II LoBALs) are even redder than quasars with absorption in high-ionization species (e.g., C IV and N V; HiBALs) on average (Weymann et al. 1991;Brotherton et al. 2001;Reichard et al. 2003). This phenomenon is readily interpreted as a result of dust extinction (e.g. Sprayberry & Foltz 1992;Voit et al. 1993;Reichard et al. 2003;Hewett & Foltz 2003;Dai et al. 2008;Jiang et al. 2013). The BAL outflow launching region is suggested to be co-spatial with or outside of the BEL regions, based on the observational fact that BAL features usually obscure BELs. Dust can survive at the outer boundary of BEL regions according to the reverberation mapping results for local active galactic nuclei (Suganuma et al. 2006). It implies the possibility of dusty outflows. Recently, Grupe et al. (2013) found an interesting anti-correlation between luminosity and UV continuum slope in a variable BAL Seyfert 1, WPVS 007. It also favors a dusty outflow component moving transversely across our line of sight. The composite spectra studies suggested an SMC-like reddening law for dust associated with quasar outflows (Reichard et al. 2003;Zhang et al. 2010), though it may not always be the case (see the studies of some individual objects, e.g. Hall et al. 2002;Jiang et al. 2013). The average amount of reddening in HiBAL and LoBAL quasars are E(B − V ) ∼0.023 and E(B − V ) ∼0.077, respectively (Reichard et al. 2003;Zhang et al. 2010). Such large amount of radiation absorbed/scattered by dust grains must be partly re-radiated in infrared band, yielding a prominent and detectable feature. More recently, Wang et al. (2013, hereafter Paper I) found an unexpected correlation between the amount of hot dust emission relative to accretion disk emission and the blueshift of C IV BEL in z ∼ 2 non-BAL quasars. The correlation dramatically strengthens with increasing Eddington ratio, since outflows tend to be dominant in high Eddington ratio quasars. It strongly implies an important role of dust in the outflow physics, such acceleration mechanisms or interaction with nearby medium (Paper I), and even in the galaxy physics and cosmology (see Elvis et al. 2002).
BALs and blueshifted BELs are very likely to be different appearances of the same outflow component viewed from different inclination angles (see e.g. Wang et al. 2011). The normal BELs emitted by virialized (or rotational) gases may hamper the proper measurement of outflow properties from BELs. Fortunately this effect is negligible in BAL measurements. It is thus necessary to revisit and extend the relationship between outflow properties and hot dust emission using a BAL quasar sample. In this paper, we construct a large z ∼ 2 BAL quasar sample from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) and the Sloan Digital Sky Survey (SDSS; York et al.2000). Following Paper I, we adopt the rest-frame NIR slope, measured from WISE data, as an indicator of the relative amount of hot dust emission. We use the BAL parameters measured from SDSS spectra, such as velocities and absorption strength, to indicate the outflow strength as commonly adopted in the literature. This paper is organized as follows. The sample construction and outflow parameter measurement are shown in Section 2. We analyze the data and present the correlations in Section 3. Finally, we discuss the possible underlying physics in Section 4, and summarize the results in Section 5. Throughout this paper, we adopt the CDM 'concordance' cosmology with H0 = 70 km s −1 Mpc −1 , Ω m = 0.3, and Ω Λ = 0.7.

SAMPLE CONSTRUCTION
We start from the SDSS DR7 quasar catalog (Schneider et al. 2010), and select 23,059 quasars with redshifts of 1.7 < z < 2.2. The redshift range is chosen to ensure reliable measurements of C IV BALs and Mg II BELs, the latter can be used to derive quasars' fundamental parameters. We coordinate cross-match these quasars with entries from the WISE All Sky Data Release catalog, using a match radius of 3 arcseconds, and get 21,720 objects detected by WISE. WISE surveyed the whole sky at four infrared bands, i.e. W1, W2, W3 and W4. The effective wavelengths of the first three bands are λ 1 = 3.35, λ 2 = 4.60 and λ 3 = 11.56, which covers the rest-frame NIR band (from 1µm to 4µm) for the selected quasars. Since we are interested in the rest-frame NIR emission, we only consider the first three bands in this work. After discarding the quasars with signal to noise ratio (S/N) less than three in the three WISE bands, we obtained 18,451 quasars that compose the DR7-WISE quasar catalog. Then cross-correlating the DR7-WISE quasar catalog with the SDSS DR5 learning vector quantization (LVQ) C IV BAL quasar catalog (Scaringi et al. 2009, hereafter SC09) and the SDSS DR7 BAL quasar catalog we compile a BAL quasar candidate catalog, which contains 2964 BAL quasar candidates (these candidates are identified as BAL quasars either by SC09 or SH11).

Identification and Measurement of Absorption
Troughs This candidate catalog is not our final working sample. The BAL quasar candidates are collected from two different samples in the literature, so the selection criterion in our candidate catalog is not self-consistent. For instance, a considerable amount of candidates are included by one work, but rejected by the other. In addition, only part of candidates have BAL parameter measurements. Therefore it is essential to use the same criterion to refine the sample, and then measure the BAL parameters for the whole sample. We develop a new method to construct our working sample from the candidate catalog. The details of our algorithm are presented in Appendix A, and here we only briefly describe it. The basic idea, same as previous works, is to construct an unabsorbed model spectrum (hereafter model spectrum) which adequately resembles the unabsorbed spectral regions of a given candidate, and then compare it with the given spectrum to identify potential absorption troughs.
For a given BAL quasar candidate, the model spectrum is obtained through the following procedure. 1) Construct a quasar template spectrum library by randomly selecting 400 quasars from a sample of SDSS non-BAL quasars with redshifts similar to the BAL quasar candidate. 2) Use each template spectrum multiplied by SMC extinction law with a free parameter E(B-V) to model the unabsorbed spectral regions of the given candidate spectrum. The parameter E(B-V) is determined by minimizing the χ 2 . 3) Rank the 400 model spectra in the increasing order of χ 2 and select the model spectrum, which has the minimum difference in the BEL regions from the candidate's spectrum, among the first 20 model spectra as our final model spectrum.
Our method can choose an model spectrum which well matches the unabsorbed continuum and BELs of the BAL quasar candidate (Figure 1). It is worthwhile noting that the SMC curve is used to adjust the template spectrum to match the given spectrum, so the E(B-V) parameter can be either positive or negative, dependent on the continuum slopes of the two spectra.
The spectra of the BAL quasar candidates are then normalized by the best model spectra. In the literature, various selection criteria were adopted to identify BALs in the normalized spectra. We use the modified "balnicity index" (BI; Weymann et al. 1991) for C IV BALs and modified "absorption index"(AI; Zhang et al. 2010) for Mg II BALs. Only C IV absorption which dips below 10% of the model spectrum and is at least 2000 km/s broad located between 0 and 30,000 km/s blueward of the quasar redshift is included in the BI calculation, while Mg II absorption at least 1600 km/s broad located between 0 and 20,000 km/s is included in AI calculation. We use different selection criteria to identify C IV and Mg II BALs, because it is well known that Mg II BALs are generally weaker and narrower than C IV BALs Trump et al. 2006;Gibson et al. 2009, hereafter G09;Allen et al. 2011, hereafter A11). In the BI calculation, the integration starts from zero velocity (see -Observed spectra (black) overplotted with the corresponding best-fitting model spectra (red) for four BAL quasars (marked with "B") and two non-BAL quasars (marked with "N"). Note that the two non-BAL quasars had been identified as BAL quasars by other works (see text for details). The gray regions are masked out in the second loop of step 2 and the calculation of χ 2 and χ 2 EL . They are considered to be the potential broad absorption regions or narrow absorption lines (see Appendix A for details). The number in each panel is the reduced χ 2 .
also Trump et al. 2006;G09;DiPompeo et al. 2012) rather than the traditional velocity of 3000 km/s (Weymann et al. 1991). Such choice can lead to a higher completeness of low velocity BAL quasars compared to previous works (see Appendix B). Additionally, according to the comprehensive study in Zhang et al. (2010), using a condition of continuous absorption width of 1600 km/s in AI calculation leads to a good trade-off between the correctness and completeness of LoBALs.
Quasars with BI > 0 and AI = 0 are classified as Hi-BAL quasars, and quasars with both larger than zero are classified as LoBAL quasars. Finally, our working sample contains 2099 BAL quasars, 264 of which are LoBAL quasars. Accordingly, 16352 quasars with BI = 0 in our DR7-WISE quasar catalog are classified into the non-BAL quasar sample. In Figure 1, we show the spectra of six BAL quasar candidates, four of which are finally classified as BAL quasars by our algorithm. We overplot the chosen model spectra for comparison. And the spectral regions masked out in the fitting process are shown in gray (see Appendix A for details). The chosen model spectra look quite similar to the observed spectra. They illustrate that our algorithm can effectively detect BALs with various profiles and velocities, e.g. multiple troughs BAL: SDSS J011605.94+133402.3, high velocity BAL: SDSS J082841.53+451233.9 and BAL superimposed to BEL: SDSS J010727.38-105743.1, and reject quasars with only narrow absorption, e.g. SDSS J004613.47+002358.1 and SDSS J022203.23-091531.3. The quasar J035310+094418.7 has quite low signal to noise ratio, so the match is little worse.
We also measure the velocity parameters of C IV BALs, i.e., the minimum, maximum and absorptiondepth-weighted average blueshift velocities (V min , V max , and V ave ). The positive velocities indicate blueshifted absorption with respect to the quasar's systemic redshift. In order to further examine the reliability of our method and measurement, we compare our working sample with other samples in Appendix B in great detail. The comparison shows that our method is at least as good as those employed in previous works. We also estimate the uncertainties of the measured BAL parameters using 500 simulated BAL quasar spectra. The simulated spectra are constructed via multiplying the randomly selected non-BAL quasar spectra by broad absorption profiles. To mimic BALs in a realistic way, the absorption profiles are extracted from randomly selected BAL quasars. The simulation shows that the typical 1σ errors for the maximum, minimum velocities and the balnicity are 3.6%, 10.9% and 6.7% respectively. Please see Appendix C for how we exactly perform the simulation and calculate the uncertainties.

Measurements of UV and Infrared Emission
We are interested in the correlation between outflow and NIR emission. All of BAL quasars in our sample have reliable magnitude measurements at three WISE infrared bands, W 1, W 2 and W 3. The redshifts of our quasars are close to 2, so the rest-frame wavelengths of the three bands all move into the NIR band, covering the wavelength region from 1µm to 4µm. We convert the magnitudes to monochromatic luminosities (L λ , where λ = W 1, W 2 and W 3), and then use a power law (L λ ∝ λ βNIR ) to fit the NIR SED in the rest-frame. The NIR spectral slope, β NIR , is thus obtained. As discussed in Paper I, hot dust emission begins to become dominant at wavelength larger than 1µm, compared to the big blue bump (Elvis et al. 1994), and host galaxy contribution is negligible for z ∼ 2 quasars (Hao et al. 2014). Therefore, the NIR spectral slope is a good indicator of the amount of hot dust emission relative to the amount of emission from accretion disk, in the sense that β NIR increases with increasing hot dust emission. In Paper I, for comparison, we also introduced an infrared to UV luminosity ratio, CF3= log(L W3 /L 2500 ), to indicate the hot dust emission relative to the disk emission. We found that β NIR shows stronger correlations with outflow properties than CF3, and suggested that β NIR is less affected by variation, inclination effect and dust extinction (see the details in Section 3 of Paper I).
The other interested parameters are the UV continuum shape and quasar fundamental parameters, such as black hole mass and Eddington ratio. Quasar UV continua are often approximated by power laws. Here, we use the power law, λF λ ∝ λ βUV , to fit the continua in two windows, [1690, 1700]Å and [2225, 2250]Å. The index β UV can be used to characterize the continuum shape. The two continuum windows are chosen to avoid strong BELs, Fe II multiplets and BALs. We use the empirical formula, presented in Wang et al. (2009a), to compute black hole masses (M BH ). The full width at half maximum (FWHM) of Mg II BELs and the monochromatic luminosity at 3000Å (L 3000 ), used for M BH estimation, are taken from SH11. Eddington ratio is evaluated as l E = L BOL /L Edd , where L Edd is Eddington luminosity and the bolometric luminosity L BOL = 5.9 × L 3000 (McLure & Dunlop 2004). Since the Mg II line profiles of LoBAL quasars can not be well determined, we estimate these parameters only for HiBAL quasars. We also estimate β NIR , β UV , M BH and l E for our non-BAL quasar sample.

Properties of BAL quasars
In Figure 2, we show the probability distributions of the measured parameters for both BAL and non-BAL quasars in our DR7-WISE quasar catalog, if available. The blue, black and red histograms represent non-BAL quasars, BAL quasars and LoBAL quasars, respectively. Generally, there are systematic differences in the shown properties between BAL and non-BAL quasars. BAL quasars, on average, have higher luminosities and Eddington ratios than non-BAL quasars, consistent with previous works (e.g., Reichard et al. 2003;Trump et al. 2006;Gangly et al. 2007; G09; A11). These results favor a radiatively driven outflow scenario. Meanwhile, an opposite and weaker trend is found in the black hole mass parameter. Because the SDSS spectroscopical survey is magnitude limited, low l E quasars tend to have high M BH (see e.g. Dong et al. 2009). Thus the dependence on M BH shown in this work might be not intrinsic but a combined effect of the dependence on l E and the false correlation between M BH and l E . However Gangly et al. (2007) found a different mass dependence. We note that the mass dependence is not detected in the SDSS DR5 low-z LoBAL quasar sample (Zhang et al. 2010).
The BI parameters are estimated in the same way as BI 0 used in G09, with the only difference being in integration range. We calculate BI in the velocity range from 0 km/s to 30,000 km/s, while G09 adopted a smaller range from 0 km/s to 25,000 km/s. Since outflows with velocities larger than 25,000 km/s are relatively rare, this difference would not significantly impact the parameter measurements as long as we have the same model spectra as theirs. It is thus useful to compare our parameters with theirs. Our BI distribution peaks at about 3000 km/s and G09's distribution peaks at the velocity range from 2000 to 4000 km/s. Meanwhile, our V min distribution peaks at about 2000 km/s, also consistent with G09. However, the V max distribution looks different from G09's. Our distribution has an apparent peak at ∼9000 km/s, while G09 presented an almost uniform distribution in the velocity range from 3000 km/s to 25,000 km/s. The difference may be partly due to the fact that our method can more effectively detect the low-velocity BALs. In Appendix B, we present a much more detailed comparison of our measurements and the other measurements available in the literature (G09 and A11).
Consistent with previous works (Weymann et al. 1991;Reichard et al. 2003;Trump et al. 2006  Correlation coefficients between β NIR and other quasar physical parameters Note. For each entry, we list the Spearman rank correlation coefficient (ρr) and the probability of the null hypothesis that the correlation is not present (P null ).
BAL quasars are, on average, redder in UV spectral region than non-BAL quasars. It suggests the existence of a significant amount of obscuring dust associated with BAL (especially LoBAL) outflows. Besides the UV continuum slope, LoBAL quasars have different BI distribution from the entire BAL quasar population. LoBAL quasars are more frequently detected in BAL quasars with stronger C IV BALs (see also G09 and Baskin et al. 2013). The mean V max of C IV BALs for LoBAL quasars is 13,136 km/s, higher than that for the entire sample, 11,785 km/s. And the mean V min for LoBAL quasars is 2903 km/s, lower than that the mean value of the entire sample, 3927 km/s (see A11 for a similar result). The differences in velocities and absorption strength imply that LoBAL outflows are stronger and of larger column density than HiBAL outflows. The distributions of the NIR slopes, β NIR , appear to be similar between non-BAL and BAL quasars (the bottom right panel of Figure 2). The mean β NIR for these two samples are 0.43 and 0.49, and the standard deviations are 0.24 and 0.23. Hot dust emission in BAL quasars is only slightly stronger than that in non-BAL quasars. It is consistent with the prediction of the BAL unifying model, in which the only difference between BAL and non-BAL quasars is viewing angle (e.g. Weymann et al. 1991). The mean β NIR for LoBAL is 0.45, similar to the other two samples, while the standard deviation is apparently larger, about 0.32. If BAL outflows are really observable along some specific directions, such as equatorial direction as suggested by previous works, the similar distributions further imply that the dependence of β NIR on inclination angle is quite weak, even absent. More recently, Runnoe et al. (2013) investigated the inclination dependence of quasar SED in a small radio-loud sample. Using the radio core dominance to indicate the orientation, they show in their figure 6 that the NIR (1 ∼ 4µm) slope of the face-on quasars is slightly less than that of the edge-on quasars, broadly consistent with the expectation from our result.

Correlations of outflow with quasar physical
properties Lots of observational works have been devoted to understanding the launch and acceleration processes of quasar outflows. The correlations of outflow properties with Eddington ratio, luminosity and SED shape have been reported in the literature (Laor & Brandt 2002;Gangly et al. 2007;Fan et al. 2009;Zhang et al. 2010;Baskin et al. 2013). These findings are largely confirmed by our new BAL quasar sample. In Table 1, we list the Spearman correlation coefficients and the probabilities of the null hypothesis for correlations between the BAL parameters (V max , V min , V ave and BI) and quasar properties (l E , L 3000 , and β UV ). One can find moderate or weak (but significant) correlations between these two sets of parameters. A quasar with higher l E or higher L 3000 or bluer UV spectrum has a tendency to harbor a stronger outflow. Generally, the velocity parameters are more strongly correlated with the quasar properties than BI. Moreover, the correlation with l E is, on average, stronger than the correlations with L 3000 and β UV . And M BH is almost uncorrelated with the velocity pa- rameters and weakly related to BI. All of these results are broadly consistent with the findings in Ganguly et al. (2007), although our correlation coefficients are slightly smaller. For further comparison, we show V max against l E , L 3000 and β UV in Figure 3. Following Ganguly et al., we also plot the predicted scaling relation between V max and l E (the normalization is taken from Ganguly's figure 6) and the V max -L 3000 fitting curve (see their figure 7 and equation 5) in the corresponding panels. These two curves appear to be the upper envelopes of these scatter points, again consistent with Ganguly et al..
Now we explore the role of β NIR in the nature of outflows. Figure 4 shows the four BAL parameters against the NIR spectral slope, β NIR . HiBAL and LoBAL quasars are shown in black and red dots respectively. We divide the whole BAL quasar sample into five equal-sized subsamples according to β NIR , then calculate the median values of β NIR and BAL parameters. The results are shown in filled green square (Figure 4). The error bars indicate the standard deviations of the corresponding parameters. Despite of large scatter, clear trends between these BAL parameters and β NIR are present. One can find that a quasar with stronger hot dust emission relative to disk emission (i.e. larger β NIR ) has a tendency to harbor a stronger outflow, which has higher maximum, minimum and average velocities and is able to yield stronger absorption. As discussed earlier, we found a strong correlation between β NIR and blueshift and asymmetry index (BAI) of C IV BELs in non-BAL quasars, especially in high Eddington ratio quasars (see paper I). If the blueshifted C IV BELs are indeed emitted by the same outflow components as those producing C IV BALs (Wang et al. 2011), one would readily expect the correlations shown in this work. Our results thus confirm the connection between outflow strength and hot dust emission from a very different aspect.
We use the Spearman rank correlation tests to evaluate the correlation strength (see Figure 4 and Table 1). As expected, these correlations are moderate and significant. We also list the correlation coefficients for Hi-BAL and LoBAL quasar subsamples separately in Table  1. One interesting phenomenon is that the correlations for LoBAL quasars are stronger than those for HiBAL quasars. It hints that dust indeed plays an unusual role in these correlations. Although the physical meaning of BAL parameters is more specific than BAI (the outflow indicator measured from C IV BELs), BAL parameters are more weakly correlated with β NIR . Note that the coefficient for the BAI-β NIR relationship is 0.55 for high Eddington ratio quasars (paper I). One possible reason is that BAL parameters are absorption properties while BAI and β NIR are both emission properties. Absorption is produced by the gas along the line of sight (LOS) so the absorption properties are only dependent on the LOS structures. In contrast, emission comes from the entire emitting gas/dust, so the emission properties are able to reflect the overall properties. The different origins in the absorption and emission properties may cause large scatter and weaken the observed correlations. The second possible reason is that outflow strengths might be affected by various physical processes. It is already known that outflow strength is significantly correlated with other physical parameters as shown in previous works and this work, and β NIR is almost independent of these factors (see below). In addition, the inclination dependence (see Section 4.1) and the uncertainty in BAL parameter measurements may also introduce scatter into the correlations .
Are the correlations between β NIR and BAL properties intrinsic or a secondary effect induced by other relationships? For instance, if β NIR tightly correlates with l E , the correlations of BAL properties with β NIR could be interpreted as being induced by the correlations with l E . It is therefore necessary to check whether β NIR correlates with other quasar physical properties that are found to be related to BAL properties, such as l E , L 3000 and β UV . We list the corresponding correlation coefficients in Table 1. As one can see that the correlation between β NIR and l E is absent, consistent with that in non-BAL quasars (paper I). And β NIR only moderately or weakly correlates with β NIR and L 3000 . Therefore the secondary effect hypothesis is unlikely. More importantly, in most cases, the correlations between BAL properties and β NIR are almost the strongest among all the correlations that we investigate (see Table 1). It also favors that the correlations of β NIR with outflow properties is intrinsic and can not be induced by the other correlations shown in this paper.
The anti-correlation between β NIR and β UV itself is interesting. In a dust extinction and re-emission scenario, a reasonable expectation is a positive correlation, which however is opposite to what we find. It means that the observed correlation is possibly ascribed to another effect. For a given UV luminosity, the disk emission of a blue quasar is weaker at NIR band than that of a red quasar. Since β NIR reflects the ratio of hot dust emission to disk emission at NIR band, β NIR tend to be larger in the blue quasar. The observed correlation is naturally yielded and does not conflict with the dust extinction and re-emission scenario. The negative trend therefore suggests that dust extinction only has little impact on the UV SED, on average. In fact, the difference in β UV distributions between non-BAL and BAL quasars is much smaller than the distribution range of β UV (the right panel of the second row in Figure 2). It is consistent with previous result that the average amount of reddening in HiBAL quasars is only E(B −V ) ∼0.023 (Reichard et al. 2003). All of these suggest that the UV slope is not determined by dust extinction but other factors, such as accretion process. We further check the trend and find that the correlation is largely dominated by a small fraction of quasars with small β NIR . After excluding quasars with β NIR < 0.1 (about 6% of BAL quasars), the correlation coefficient for β NIR -β UV relationship is reduced to about 0.1.

Composite spectra
As shown in the above subsection, we present the correlations of outflow properties with quasar physical properties. There is large scatter in these correlations and the correlation strength is only moderate or weak (however significant). One possible reason for this is that the outflow properties measured from the absorption troughs are sensitive to the LOS structures. If outflows are clumpy as suggested by previous works (e.g. Smith & Raine 1988;Netzer 1993;Krolik & Kriss 1995;Netzer 1996;Elvis 2000;Chelouche & Netzer 2005), the BAL parameters may significantly vary with the viewing angle. In other works, for a given quasar outflow, the BAL parameters may be very different along different LOS. One useful method to suppress this LOS effect is to construct composite spectra. Composite spectra may eliminate the variation in outflow properties and distinctly tell us how the mean absorption profile at different velocity changes with different physical parameters.
The whole BAL quasar sample is divided into five equal-sized subsamples, according to the UV slope, Eddington ratio (here, only HiBAL quasars are used), the luminosity at 3000Å and the NIR slope, respectively. We adopt the method of Vanden Berk et al. (2001) to create the geometric mean composite spectra for each subsample. For each BAL quasar, we deredshift its observed spectrum into its rest frame. The spectrum is then normalized at 1700Å, rebinned into the same wavelength grids, and geometrically averaged bin by bin. In the left panels of Figure 5, we plot composite spectra of these subsamples in sequence. The SDSS DR7 composite spectrum is overplotted in gray dashed line for comparison. We then use the exactly same method as shown in Section 2.1 and Appendix A to obtain the model spectrum for each composite spectrum. The normalized composite spectra around C IV BALs are shown in the right panels of Figure 5 and the measured BAL parameters (V max , V min , V ave and BI) as functions of β UV , l E , L 3000 and β NIR are shown in Figure 6.
One can see, from Figure 5, that the dependence of the mean absorption profile on β NIR is the strongest. The absorption shifts blueward as a whole with increasing β NIR . Both the maximum and minimum velocities significantly increase with the increase of β NIR . And the absorption profile apparently broadens as β NIR increases. We also check Si IV BALs and find the same trend. In particular, the trend for Si IV is almost as strong as that for C IV. The second most important factor is Eddington ratio. Eddington ratio seems to have a dramatic impact on the large velocity part of outflows. The maximum velocity increases with increasing Eddington ratio. Interestingly, the absorption profile almost remains unchanged at velocity less than 10,000 km/s as l E varies. On the contrary, β UV only affects the low velocity part of the outflowing gases. The absorption trough in a red quasar is closer to zero velocity than that in a blue one. And the absorption profile is almost independent of L 3000 . The -Composite spectra of BAL quasar subsampes according to Eddington ratio, the UV slope, the luminosity at 3000Å and the NIR slope, and their normalized flux of C IV BALs. Please see the text on how we construct the composites. For comparison, the SDSS DR7 composite spectrum is shown by gray dashed curve in left panels, we also show the lines of the velocity V = 0 km/s, the normalized flux f N = 1 and f N = 0.9 in right panels.
dependencies of the mean absorption profile on quasar physical parameters may offer valuable insight into the launch and acceleration of outflows.
The impression obtained from the composites is confirmed by the quantitative results shown in Figure 6. Outflow velocities and strength measured from the composites increase with the enhancement of black hole accretion and hot dust emission. In particular, β NIR is the factor which has the strongest impact on BI. The parameter BI for the largest β NIR composite is more than twice as large as that for the smallest β NIR composite. The maximum velocity is more sensitive to l E than other factors. V max increases from ∼12,000 km/s to ∼20,000 km/s. And β NIR and l E appear to be both dominant for V ave . In the case of V min , β UV becomes as important as β NIR . One can see that V min changes from 0 km/s to nearly 900 km/s with decreasing β UV and increasing β NIR . Overall, the dependence of outflow properties on β NIR is more significant. We note that these results are well consistent with the statistical results shown in the above subsection (see Table 1).
Recently, Baskin et al. (2013) also analyzed the av-erage absorption profiles of DR7 BAL quasars. From their figure 9 and 11, one can clearly see that the UV slope (measured in the 1700-3000Å range) only impacts the small velocity part while the Eddington ratio primarily controls the large velocity part. Additionally, they found no significant dependence on the luminosity. These trends are well consistent with what we show here. However their results also have apparent difference from ours, in particular in the behavior associated with the UV slope. They found that the low velocity part of the absorption deepens as the UV continuum becomes redder. We do not detect this trend in our composite. The difference may be due to several reasons. First, we use very different method to construct the composite. Second, we use very different method to derive the normalized spectra. And third, about 17% BAL quasars in SH11 sample (used by Baskin et al.) are excluded by our algorithm and additional 117 quasars are added in(see Table A1). Another interesting result found by Baskin et al. is the significant dependence of absorption profile on the equivalent width (EW) of He II BEL, which is interpreted in terms of the dependence on ionization SED. It will be interesting to check whether there is any correlation between EW(He II) and β NIR . He II BELs are blended with O III] and strong C IV BELs, so it is hard to properly decompose He II BEL. We therefore use the integrated flux between 1620Å and 1650Å to calculate EW(He II), which is used in Baskin et al. (Please see their description of the method for details). The EW(He II) measured from the composites of different β NIR are 5.15, 5.56, 5.35, 5.91 and 5.55Å, in sequence of increasing β NIR . The NIR slope is apparently uncorrelated with EW(He II). To further verify the result, we also measure EW(He II) for all individual BAL quasars. The Spearman correlation coefficient for these two parameters is only -0.01 (the null probability is 0.78).

DISCUSSION
In Paper I, we adopted the blueshift and asymmetry index (BAI) of C IV BELs to indicate the outflow property and velocity, and found an interesting correlation between BAI and β NIR . In this paper, the blueshifted velocities and absorption strength of BALs instead of BAI are used to characterize the outflow strength. Again, we find significant correlations with β NIR , well consistent with the results shown in Paper I. Using outflow parameters measured from BELs (e.g. BAI) and BALs (e.g. V max , V ave and BI) both have advantages and shortcomings. On one hand, BELs are the integral of emission over entire volume of the outflow and represents the overall properties, while BAL troughs only hold the information of outflow along the LOS and thus are sensitive to clumpy structures. It may be the reason for that β NIR is more strongly correlated with BEL parameters than BAL parameters. On the other hand, the physical meanings of BAL parameters are more straightforward and specific than those of BEL parameters. The emission from virialized gas also contributes to BELs so that it is generally hard to reliably extract outflow parameters from BELs. In contrast, the contamination in BALs from other line features is usually unimportant and easily to handle (see Section 2). It is one of the reasons that we think it is essential to use BAL quasars to revisit the relationship between outflow and hot dust. Overall, these two works together present robust evidences for the important correlation and connection between outflow strength and hot dust emission.
In order to understand the origin of these correlations, we present detailed discussions in the following. In Section 4.1, we examine whether or not a dust-free outflow, together with a hydrostatic torus model, can properly accommodate our findings. Our analysis based on current observational results suggest that this scenario is likely unviable. In Section 4.2, we discuss the dusty outflow scenarios, which have ever been proposed in Paper I.

Dust-free outflow scenario
Suppose that hot dust emission in these high-redshift quasars has nothing to do with outflows, and is predominantly emitted by the innermost part of a hydrostatic and optically thick torus, as suggested by lots of previous studies (e.g., Neugebauer et al. 1987;Barvainis 1987;Suganuma et al. 2006;Kishimoto et al. 2007;Mor et al. 2009;Mor & Trakhtenbrot 2011). It would be interesting to see whether or not the correlations shown in this paper and Paper I are induced by a third factor that simultaneously governs or relates to outflows and dust emission.
We first discuss the inclination effect. The NIR slope, β NIR , measures the amount of hot dust emission relative to accretion disk emission. The disk emission is theoretically predicted to scale with the cosine of the disk inclination angle (IA). If hot dust emission is (approximately) isotropic, β NIR is expected to increase with increasing IA. However, dust torus may block the hot dust emission at a large IA (Roseboom et al. 2013). It complicates the situation and makes the dependence of β NIR on IA uncertain. More recently, Runnoe et al. (2013) used a sample of radio-loud quasars to investigate the inclination dependence of quasar spectral energy distribution. The composites in their figure 6 clearly present a weak dependence of β NIR on IA, in the sense that a edge-on quasar spectrum tend to show a slightly larger β NIR than a face-on one. We thus still assume that β NIR increases with increasing IA Since the outflow covering factor is usually small, about 10%∼20% (e.g., Tolea et al. 2002;Hewett & Foltz 2003;Reichard et al. 2003;Trump et al. 2006), it is unlikely that outflow properties and β NIR systematically and significantly vary with IA within the opening angle of outflows. However, one can still elaborately design a outflow model to match the observation. For example, the mean outflow IAs vary among quasars and the small-IA outflows are, on average, weaker than the large-IA outflows. The problem for this hypothesis is that there is no any observational evidence for it and the inclination dependence of β NIR is actually weak (Runnoe et al. 2013). In addition, it also fails to account for the correlation with BAI, the outflow parameter measured from BELs (see Paper I). The inclination effect might be the other factor for the large scatter in the observed correlations.
Lots of studies have found the correlations of outflow properties with quasar fundamental parameters e.g., Reichard et al. 2003;Gangly et al. 2007;Fan et al. 2009;A11;Baskin et al. 2013). Outflow strength is significantly dependent on Eddington ratio, broad band SED (e.g. ionization SED, and UV continuum slope, β UV ) and luminosity. If these factors also have an impact on hot dust emission, one might find a piece of evidence to support the dust-free outflow scenario. As shown in this paper and previous studies (e.g. Mor & Trakhtenbrot 2011; Paper I), both β NIR and the infrared to UV flux ratio are almost independent of Eddington ratio. And the correlations between luminosity and the two hot dust indicators are either absent or negative. In addition, the correlation of outflow strength with β NIR is stronger than that with β UV and the correlation between the two slopes are apparently weaker (Table 1). Hence none of them is the factor that we are seeking. The role of ionization SED is unclear, since we do not know its relation with β NIR . Recently, Baskin et al. (2013) detected an interesting dependence of outflow strength on EW(He II) and ascribed it to the impact of ionization SED. We thus check the relationship between EW(He II) and β NIR . Again, no correlation is found.
Another interesting factor is metallicity. Quasars harboring strong outflows tend to have high gas metallicity (Hamann 1998;Leighly 2004;Wang et al. 2009b). Recently, Wang et al. (2012) found that C IV blueshift increases with gas metallicity. It means that outflows are stronger in higher metallicity environment. Since dust forms more easily in higher metallicity gas, one may expect that metallicity is an appropriate factor which simultaneously influences outflow and hot dust emission. However, such naive expectation may not be true. For a typical torus, the relative amount of dust emission is determined by dust covering factor, but not dust amount. Theoretically, it is unclear how to construct a relationship between dust covering factor and metallic-ity. And observationally, there is no evidence supporting it. On the contrary, the required positive correlation seems inconsistent with the current observational facts that metallicity strongly increases with increasing luminosity(e.g. Hamann & Ferland 1999;Nagao et al. 2006) and dust covering factor (indicated by the infrared to UV flux ratio) is weakly anti-correlated with luminosity (e.g. Mor & Trakhtenbrot 2011).
To sum up, the current observational results listed in this paper and the literature do not favor the dust-free outflow scenario. However further works, especially on the ionization SED and metallicity, are still required to examine our point.

Dusty outflow scenarios
In Paper I, we attempted to propose two plausible mechanisms to interpret the correlation of outflow strength with hot dust emission. Here we present further discussion. In the first mechanism, dust is intrinsic to outflows. This idea is supported by the similar locations of BAL outflows and hot dust (Elvis 2000 and references therein). There is little observational constraint on the (relative) location of BAL outflows. Since the absorption material almost always absorbs the BELs, it is usually believed that BAL outflows are co-spatial with or outside of BEL regions (BLRs). The reverberation mapping results of local AGNs suggested that hot dust is also located at the boundary of BLRs (e.g., Suganuma et al. 2006). These two implies the possibility of originally dusty outflows. Outflows may emerge from the outer region of accretion disk or even the innermost region of torus, in which the gas clouds are dusty and relatively cold. These dusty clouds are uplifted above the disk, and exposed to the central engine. The low density part is highly ionized and responsible for the blueshifted absorption and emission lines. Dust survives in the dense region and radiates in NIR band. In addition, dust perhaps forms in dense clouds embedded in outflows (Elvis et al. 2002). A strong outflow carries large amount of dust, and thus enhances the NIR emission. At the same time, the acceleration due to dust absorption and scattering make the outflow stronger, which further enhances the correlation.
It is worth noting that Elitzur & Shlosman (2006; see also Emmering et al. 1992) presented a torus outflow model. In their model, dusty molecular clouds are uplifted from the accretion disk and injected into outflows, and these moving optically thick clouds constitute a torus. Our scenario resembles theirs. But there are apparent differences. Firstly, the velocity of the outflows we studied is of order 10000 km/s, much faster than the torus outflows. Second, the clouds in BAL outflows are unlikely to be optically thick in dust opacity, since a large fraction of BAL quasars can be detected at UV band and the average amount of reddening in BAL quasars is only E(B − V ) ∼ 0.023 (Reichard et al. 2003). Thirdly, the large fraction of BAL outflows is highly ionized. One solution to the discrepancies is that the torus structure is stratified. The dusty BAL outflow is perhaps the inner edge of the torus outflow of Elitzur & Shlosman and is in an extreme ionization and dynamical state compared to the bulk torus. In this case, it is reasonable to speculate the existence of an outflow component with ionization and velocity in between. Recently, Zhang et al. (2013) found, in low redshift AGNs, an interesting connection between mid-infrared (MIR) emission and outflows in NLR, the velocity of which is of order several hundred kilometers per second (see also Hönig et al. 2013). The NLR outflows seem to fill the gap. However, a simple extrapolation from local AGNs may not be appropriate. Studies of the relationship between various outflow components (such as NLR, BAL, BEL) and MIR emission of z ∼ 2 quasars may help in uncovering the underlying physics.
Before discussing the second mechanism, let us list several observational facts. The maximum velocity of a typical BAL outflow is about 10000 km/s. The typical width of BELs is about 4000 km/s and BAL outflows are cospatial with or outside of BEL regions. It means that the typical Keplerian velocity at BAL outflow location is less than 4000 km/s. The last fact is the absence of evidence for outflows with similar velocity at galactic scale or in NLR. All of these facts together show that BAL outflows are strongly decelerated before reaching NLR and the gravity of their central black holes, however, is not responsible for the deceleration. Thus it is likely that outflows interact with surrounding medium. This medium is located between BLR and NLR and should be dense enough that it can effectively decelerate the outflows. Apparently, interaction with torus clouds is the most favorable.
Recent hydrodynamical simulation (Wagner et al. 2013) showed that fast outflows can break dense clouds into diffuse warm filaments. This process makes more dust in the clouds exposed to the central UV source. As a consequence, the infrared emission increases and the outflows become dusty. Since a stronger outflow can ablate dense clouds more effectively, the dependence of NIR emission on outflow strength is yielded. In order to reach a high enough temperature to radiate in NIR band, the interaction has to occur at the innermost region of torus. In this scenario, outflows confine the geometry and subtending angle of dusty torus. It is in line with the suggestion that outflows blow away ambient dusty gas and transform a buried quasar into normal quasar phase (e.g. Sanders et al. 1988;Rupke & Veilleux 2013). One problem of this scenario is that the interaction timescale (see e.g. Wagner et al. 2013) is much shorter than the quasar lifetime. The correlation might disappear after outflows blow away all of clouds in the outflow direction.

SUMMARY
Based on a large z∼2 BAL quasar sample built from the Wide-field Infrared Survey and the Sloan Digital Sky Survey data Release Seven, we investigate the correlations of outflow strength on hot dust emission (indicated by the near-infrared continuum slope, β NIR ) and other quasar properties, such as Eddington ratio, luminosity and UV continuum slope. We collect BAL quasars from two existed BAL quasar samples in the literature, then we use our new automatic algorithm to refine the sample and measure the BAL parameters. The final working sample contains 2099 BAL quasars.
We first reexamine the differences between BAL and non-BAL quasars and the dependencies of outflow properties on Eddington ratio, luminosity and UV slopes, which have been reported in the literature. BAL quasars, on average, have redder UV continua, higher Eddington ratios and luminosities than non-BAL quasars. LoBAL quasars have more severe and slightly broader C IV absorption than HiBAL quasars. Outflow velocities and strength are found to be moderately dependent on Eddington ratio and weakly dependent on luminosity and UV slope. All of these results are consistent with previous works.
We find moderate and significant dependencies of outflow velocities and strength on β NIR , in the sense that outflows strengthen with increasing hot dust emission. It is consistent with the correlation of the blueshift of C IV BELs with β NIR in non-BAL quasars (see paper I). The large scatter in the relationship may be partly due to the fact that the outflow velocities and strength measured from absorptions are sensitive to the LOS structure. Our statistical analysis and composite spectra study both reveal that outflow strength is more strongly correlated with β NIR than Eddington ratio, luminosity and UV slope. In particular, the composites reveal that the entire C IV absorption profile shifts blueward and broadens as β NIR increases, while Eddington ratio and UV slope only affect the high and low velocity part of outflows, respectively.
We discuss several potential processes that are possibly responsible for the observed correlations of outflow properties on hot dust emission. In a dust-free outflow and hydrostatic and optically thick torus scenario, we attempt to use inclination effect and other observed relationships associated with outflows to interpret the correlations. However, our analysis does not favor this scenario. We thus suggest that dust is intrinsic to outflows and may have a nontrivial contribution to the outflow acceleration.
State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton Uni-versity, the United States Naval Observatory, and the University of Washington. adequately resembles that of a given BAL quasar. The non-BAL quasar spectrum can be used as an unabsorbed model spectrum (hereafter model spectrum) to characterize the intrinsic continuum and BELs of the BAL quasar.
The absorption troughs, if present, are then shown in the spectrum of the BAL quasar normalized by the model spectrum.
For a BAL candidate in consideration, we employ our new automatic algorithm to seek the model spectrum through the following three steps.
1. Construct a quasar template spectrum library: Four hundred quasar spectra in the redshift range of 0 < z − z 0 < 0.05 are randomly selected from the non-BAL quasar catalog (see Section 2) to build up the quasar template spectrum library for the given BAL quasar candidate. Here z and z 0 are the redshifts of the non-BAL quasar and the given candidate. Then these template spectra are smoothed by twice iteration B-spline analysis 4 . Pixels with flux intensities smaller than 85% of the smoothed spectra may have narrow absorption lines, and therefore are masked out from the following fit.
2. Model the given BAL candidate spectrum: The template spectra are scaled and multiplied by the SMC extinction law with a free parameter E(B-V) to fit the observed spectrum of the given BAL candidate in two-step fitting loops. The fitting process is performed using the MPFIT package (Markwardt 2009). The scale parameter and E(B-V) are determined by minimizing the χ 2 . After the first loop, we mask out the potential broad absorption regions where the observed flux are lower than 85% of the scaled and reddened template, and refit the spectrum. The overlap spectral regions between the candidate spectrum and the template spectra, except those masked, are used in the fitting process and the calculation of χ 2 .
3. Choose the best unabsorbed model spectrum: We rank the 400 model spectra in the increasing order of χ 2 (calculated in the second loop of step 2), and take out the first 20 model spectra. We calculate the emission line chi-square, χ 2 EL , in the wavelength range of [1450, 1650]Å, where the pixels are not masked. This parameter measures the similarity between the observation and the model in the unabsorbed C IV BEL region. The best model spectrum is chosen by minimizing χ 2 EL among the twenty model spectra. Therefore our best model spectrum has the minimum difference in the C IV BEL regions from the candidate's spectrum, and is also in good agreement with the observed broad band continuum.
For illustration, we show, in Figure 1, the masked regions in gray for several candidate quasars. Almost all of the significant narrow absorption lines and the potential broad absorption lines are well identified. The quasar in the left bottom panel has low signal to noise ratio and many individual absorptions, so a considerable large part of the spectrum is masked. In order to demonstrate the reliability of our algorithm, we show the probability distributions of the reduced χ 2 and χ 2 EL for the chosen best model spectra of all candidate BAL quasars in the right panel of Figure  A1. The typical reduced χ 2 is 1.5. Considering that we model the candidate spectra in a very large spectral region, our fitting results are quite good. We then show, in the left panel of Figure A1, the number count distribution of the reduced χ 2 of 400 template spectra which are used to model a typical candidate quasar SDSS J010727.38-105743.1.
One can see that about half of template spectra are able to properly model the given spectra. And the distribution of the reduced χ 2 EL for the 20 smallest-χ 2 model spectra are shown in red. Obviously, the differences between the 20 model spectra and candidate spectrum in C IV BELs are also very small.

B: RELIABILITY OF BAL IDENTIFICATION
To examine the reliability of our BAL quasar identification, we cross-match our working BAL quasar sample with the BAL quasar catalogs constructed by SC09, A11 and SH11. Note that the comparisons with SC09 and A11 are restricted to SDSS DR5 and DR6, respectively. Table A1 lists the numbers of BAL quasars identified or rejected by our work and the other three samples. There are 1606 learning vector quantization (LVQ) BAL quasars, 1769 hybird-LVQ BAL quasars, 1315 A11's BAL quasars and 2395 SH11's BAL quasars which also belong to our z ∼ 2 DR7-WISE quasar sample (marked with "T"). Most of them are classified as BAL quasars by this work (marked with "B", see Table A1). Moreover we reject 284 LVQ BAL quasars, 134 hybird-LVQ BAL quasars, 65 A11's BAL quasars  and 406 SH11's BAL quasars (marked with "L"). In our working sample, there are 315, 2, 621 and 119 objects which are not included in LVQ, hybrid-LVQ, A11's and SH11's BAL quasar samples (marked with "A"). The correctness and completeness of our working sample are very similar to the hybrid-LVQ BAL quasar sample identified by SC09.
In Figure A2, we compare the composite spectra of the BAL quasar subsamples listed in Table A1. The black, green and red lines represent the BAL quasar subsamples marked with "B", "L" and "A", respectively. One can see that there are almost no significant BALs in the composite spectra of the BAL candidates that we rejected but accepted by other works (green lines). On the contrary, the BALs are remarkable in the composites of BAL quasars we accepted but rejected by other works (red lines). A11 and SH11 tends to miss BAL quasars with low-velocity absorption.
As shown in Table A1, there are 1250 objects which are included in both our and A11's samples. The BAL parameters are available for these A11 objects. Moreover, among the 1982 BAL quasars included in both our and SH11's samples, 1517 objects have available BAL measurements from G09. Before presenting a comparison with these two works, we have to emphasize the differences in BAL identification among the three works, which may also have an impact on the BAL parameter measurements as shown below. As discussed in Section 2.1, BI or its modifications are the most widely used metric to identify BALs, and all of the modifications follow the formula presented in Weymann et al. (1991). The only difference is in integration range. A11 adopted the traditional integration range from 3000 km/s to 25000 km/s. G09 extend to low-velocity end to zero in the calculation of BI 0 . We adopt a even larger integration range from 0 km/s to 30000 km/s. The direct comparison among these three works is shown in the left and middle panels of Figure A3. Our parameter measurements are well consistent with those of A11 and G09, although there is apparent discrepancy and scatter. Part of the discrepancies are due to the difference in BI definition as discussed above. For example, the high and low-velocity cutoffs in the V max distribution of A11 sample are apparently ascribed to the chosen integration range. In addition, our BI parameter tends to be greater than A11's at low-BI end. It is also partly due to their smaller integration range than ours. The considerable scatter is also expected, since we use very different method to construct the model spectra.
In right panels, we show the parameter distributions of 1036 quasars, which are identified as BAL quasars by all the three works. As discussed in Section 3.1, the V max distributions are different between ours and G09. Compared to our distribution, theirs is much flattened. However, our V max distribution well matches A11's, except at the high and low-velocity ends. The discrepancies at two ends are again caused by the different integration ranges. The three V min distributions are quite similar at V min > 4000 km/s. Since the starting velocity of BALs is set to 3000 km/s, A11's distribution has a significant high peak at V min = 3000 km/s. The bottom right panel shows that our BI distribution is similar to G09's and A11' balnicity trends to smaller than ours because of the integration range.
Overall our comprehensive comparison suggests that our BAL quasar sample is at least as good as, or even better than, those in the literature.

C: UNCERTAINTIES OF BAL PARAMETERS
The measurement uncertainties of BAL parameters arise from the statistical noise and the uncertainties in model spectrum. The measurement uncertainty due to the statistical noise is usually easily to handle. However, the model uncertainties is hard estimated. We therefore construct 500 simulated BAL quasar spectra to quantify the uncertainties. The simulated spectra are constructed via multiplying the randomly selected non-BAL quasar spectra by absorption profiles. To mimic BALs in a realistic way, the absorption profiles are also extracted from randomly selected BAL quasars. Then we fit the simulated spectra and measure the BAL parameters following our automatic algorithm. We emphasize that the non-BAL spectra selected to create simulated spectra are rejected during the fitting process.
As a simple demonstration, we show, in Figure A4, the non-BAL quasar spectrum of SDSS J122626.22+270437.0 (gray curve in panel (a)), the observed C IV absorption profile of SDSS J035310.23+094418.7 (red line in panel (b)), and the simulated BAL quasar spectrum (black curve in panel (a)). The best model spectrum for the simulated spectrum is shown in red line in panel (a) and the normalized flux in gray line in panel (b). The identified absorption profile (green line in panel (b)) is very similar to the input profile, despite of tiny difference. One can see that the maximum velocity is reduced because of the flux noises. And the absorption close to 1549Å is slightly shallower than the input one, which is likely caused by the BEL profile mismatch.
We have original V max , V min and BI for the 500 simulated spectra. These original quantities, X input , are directly measured from the input profiles. We also get the corresponding parameters from the simulated spectra, X sim , the absorption profile identified from the simulated spectrum. Gray line: the normalized flux of the simulated spectrum. Panel (c-e): the distributions of the relative difference ∆x between the input parameters and the parameters measured from the simulated spectra. Here ∆x = (x input − x sim )/x input , x in the three panels denotes Vmax, V min and BI respectively. using our automatic algorithm. So we can derive the uncertainties in the three parameters using the formula ∆ x = (x input − x sim )/x input . The distributions of ∆ x are shown in panel (c-e) of Figure A4. The standard deviations of the ∆ x distributions can be used to represent the typical 1σ errors. The 1σ errors for V max , V min and BI are 3.6%, 10.9% and 6.7%, respectively. We note that the resultant errors contain both the uncertainties due to statistical noise and the model uncertainty.