domain: phys shs

text:
Vincent Verbavatz

Marc Barthelemy
email: marc.barthelemy@ipht.fr

Critical factors for mitigating car traffic in cities

Car traffic in urban systems has been studied intensely in past decades but its analysis is often limited to empirical observations and agent-based modelling, and despite the importance and urgency of the problem we have a poor theoretical understanding of the parameters controlling urban car use and congestion. Here, we combine economical and transport ingredients into a statistical physics approach and propose a parsimonious model that predicts the share of car drivers, the CO2 emitted by cars and the average commuting time. We confirm these analytical predictions on 25 major urban areas in the world, and our results suggest that urban density is not the relevant variable controlling car-related quantities but rather are the city's area size and the density of public transport. Mitigating the traffic (and its effect such as CO2 emissions) can then be obtained by reducing the urbanized area size or, more realistically, by improving either the public transport density or its access. In general, increasing the population at fixed area would increase the emission of CO2 in sharp contrast with the commonly accepted paradigm that increasing the density leads to a reduction of car traffic.

INTRODUCTION

As most humans now live in urban areas and two-third of the world population will live in cities by 2050 [START_REF]World Urbanization Prospects: The 2014 Revision[END_REF], understanding urban mobility patterns [START_REF] Barbosa | Human mobility: Models and applications[END_REF] has become paramount in reducing transport-related greenhouse gas emissions and crucial to efficient environmental policies [START_REF] Dodman | Blaming cities for climate change? An analysis of urban greenhouse gas emissions inventories[END_REF][START_REF] Glaeser | The greenness of cities: Carbon dioxide emissions and urban development[END_REF][START_REF] Newman | The environmental impact of cities[END_REF][START_REF] Oliveira | Large cities are less green[END_REF]. In a seminal paper, Newman and Kenworthy correlated transport-related quantities (such as gasoline consumption) with a determinant spatial criterion: urban density [START_REF] Newman | Cities and automobile dependence: An international sourcebook[END_REF]. Higher population density areas were shown to have reduced gasoline consumption per capita and thus reduced gas emissions. Their result had a significant impact on urban theories over the last decades and has become a paradigm of spatial economics [START_REF] Creutzig | Global typology of urban energy use and potentials for an urbanization mitigation wedge[END_REF]. This study is however purely empirical and has no theoretical foundation, which casts some doubts about the importance of density as the sole determinant of gasoline consumption and other car dependent quantities.

On the other hand, car traffic has been studied at various granularities and with various tools ranging from agent-based modelling (see for example [START_REF] Balmer | Towards truly agent-based traffic and mobility simulations[END_REF][START_REF] Bazzan | A review on agent-based technology for traffic and transportation[END_REF]), cellular automata [START_REF] Nagel | A cellular automaton model for freeway traffic[END_REF], or hydrodynamic approaches (see for example the review [START_REF] Nagatani | The physics of traffic jams[END_REF] on various physics type approaches). Car congestion was also considered from the point of view of economics (see for example [START_REF] Arnott | Economics of a bottleneck[END_REF][START_REF] Duranton | The fundamental law of road congestion: Evidence from US cities[END_REF]). We have however no parsimonious model that is able to point to the critical parameters and dominant mechanisms of traffic in urban areas. In this paper, we address this problem and propose a theoretical model of urban daily commuting relying on two crucial ingredients: the coexistence of cars and mass rapid transit (MRT), and traffic congestion. Combining these ingredients within a disordered system type approach allows us to derive significant conclusions concerning MRT transport ridership and transport-related greenhouse gas emissions. We compare our predictions to empirical data obtained for 25 major cities in the world (see Sup. Mat. for details about data) showing an excellent agreement given the simplicity of the model and the absence of any ajustable parameter. In this approach we deliberately left out details of these systems and focused on the basic processes that capture the complexity of urban systems while accounting for qualitative and quantitative behaviors [START_REF] Pumain | Scaling laws and urban systems[END_REF][START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF][START_REF] Bettencourt | A unified theory of urban living[END_REF][START_REF] Batty | The size, scale, and shape of cities[END_REF][START_REF] Barthelemy | The structure and dynamics of cities[END_REF].

THE MODEL

According to the classical urban ecomomics model of Fujita and Ogawa [START_REF] Fujita | Multiple equilibria and structural transition of non-monocentric urban configurations[END_REF], individuals choose job and dwelling places that maximize their net income after deduction of rent and commuting costs. More precisely, an agent will choose to live in x and work at location y such that the quantity

Z(x, y) = W (y) -C R (x) -G(x, y) (1)
is maximum. The quantity W (y) is the typical wage earned at location y, C R (x) is the rent cost at x, and G(x, y) is the generalized transportation cost to go from x to y. There is also a similar equation for the profit of companies (that they want to maximize) and that we do not need here: we assume that employment is located at a unique center y = 0 and that wages and rent costs are of the same order for all individuals. Most large cities are polycentric [START_REF] Louf | Modeling the polycentric transition of cities[END_REF][START_REF] Louf | How congestion shapes cities: From mobility patterns to scaling[END_REF] with the existence of many different activity centers but this first approach assumes that polycentricity does not change the order of magnitude of commuting trip distances (see also [START_REF] Louf | How congestion shapes cities: From mobility patterns to scaling[END_REF] on this point).

We also assume that the residence location x is given and random -residence choice is obviously a complex problem and replacing a complex quantity by a random one is a typical assumption made in the statistical physics of complex systems. Within these assumptions, we obtain a simplified Fujita-Ogawa model where the maximization of Z(x, 0) is equivalent to the minimization of the transport cost G(x, 0). Thus, individuals simply tend to minimize their commuting costs and we will discuss their commuting trips between their home and the central business district. All these assumptions are of course approximations to the reality but we claim here that our model captures the essence of the traffic in large urban areas phenomenon. Starting with a model containing all these various parameters would actually be not tractable and would hide the critical ingredients. We assume that a proportion p of the population has access to mass rapid transit such as the subway or elevated rail (we neglect buses or tramways here) whereas a share 1p of the population has no choice but to commute by car (we assume that all individuals can drive a car if needed). We omit all other forms of commuting (walking, cycling, etc.), and we neglect spatial correlations between the densities of public transport and residence or population, which is an important assumption that certainly needs to be refined in future studies. The fraction p of individuals that have a choice between car and MRT will choose the transport mode with the lowest generalized cost which takes into account both monetary costs and trip duration (see for example [START_REF] Glaeser | Why do the poor live in cities? The role of public transportation[END_REF][START_REF] Crozet | Le temps et les transports de voyageurs (CEMT-ECMT)[END_REF]). For cars, we include congestion described by the Bureau of Public Roads function (see [START_REF] Branston | Link capacity functions: A review[END_REF] and Methods) which captures the main effect: increasing the traffic on a road will decrease the effective speed on it. The corresponding generalized costs for cars and the MRT then read

G car (x) = C c + d(x) v c V 1 + T c µ (2)
G M RT (x) = V f + d(x) v m (3)
where C c is the daily cost of a car, v c and v m are respectively the car and MRT velocities, c is the road capacity of the city, f the walking plus waiting time for transit, d(x) the distance between home located at x and the central business district, T the total car traffic, and µ the exponent that characterizes the sensitivity to traffic. The quantity V is the value of time defined in transport economics as the money amount that a traveler is willing to pay in order to save one hour of time. It is an increasing function of income and is bounded by the hourly wage. We assume here that driving is faster than riding public transport (v c > v m ) but is more expensive (C c > 0we neglect here the cost of a MRT ticket in comparison with car costs). Once an individual has chosen a mode he sticks to it and will not reconsider his choice even if the traffic evolves. In other words, we assume here that individual habits have a longer time scale than traffic dynamics.

Individual mobility is then governed by comparing these costs G car and G M RT and will depend on exogenous parameters (such as car and subway velocities, car costs, etc.) and endogenous parameters such as the commuting distance. In the general mode choice theory (see for example [START_REF] Lerman | Discrete choice analysis: theory and application to travel demand[END_REF]), given the values of the costs G car and G M RT there is a probability P C = F (G car -G M RT ) to choose the car. The function F is in general smooth and satisfies F (-∞) = 1 and F (+∞) = 0 and we consider here the simplest case where F (x > 0) = 0 and F (x < 0) = 1. An individual located at x with access to the MRT will then choose to use the car if G car < G M RT which implies a condition on the value of time of the form V < V m (x) where V m (x) depends on the parameters of the system and on d(x) and reads

V m (x) = C c f + d(x) 1 vm -1 vc 1 + T c µ (4)
Also, writing the equality G car = G M RT between generalized costs of car and MRT leads to the critical distance given by

d(V, T ) = min L, Cc V -f 1 vm -1 vc (1 + ( T c ) µ ) (5)
where L ∼ √ A (A is the area size of the city) is the largest extent of the city. This critical distance evolves as traffic increases (Fig. 1), translating the fact that driving is less advantageous when the traffic is large. This expression also shows that individuals with a small value of time are more likely to use public transport, since they are more apt to spend time than money. Distance to the center is pivotal in this decision process: too far from the centerfurther than a critical distance d(V, T ) -individuals favor driving to avoid lengthy journeys, and the richer they are, the smallest this distance.

Writing d(V, T = T * ) = L, gives the critical maximal traffic T * for which driving is less beneficial in the whole agglomeration:

T * = c v c 1 v m - 1 v c - Cc V -f L 1/µ (6)
Computing the car share

In the following, we will then discuss the traffic dynamics in the two regimes (T > T * or T < T * ). For T < T * , the evolution of the car traffic when the population P increases is given by the differential equation:

dT dP = 1 -p + p 1 -π d(V, T ) 2 A (7)
where the first term of the r.h.s. corresponds to the 1p share of individuals far from MRT stations, and the second term to individuals living further from the center than d(V, T ) (we consider here statistically isotropic cities). Plugging in the expression Eq. 5 of d(V, T ) leads to an equation that can be solved exactly at least for µ = 1 and µ = 2. However, for all µ, we obtain at dominant order in

T /c T 1 - p Ab 2 P + O T c µ ( 8
)
where b is a function of the exogenous parameters

b = 1 v m - 1 v c 1 √ π(C c /V -f ) (9)
For P > P * = T * / 1 -p Ab 2 , since d(V, T ) > L the only source of car traffic comes from individuals who do not have access to the mass rapid transit, which leads to dT /dP = 1p implying T = (1p)(P -P * ) + T * [START_REF] Bazzan | A review on agent-based technology for traffic and transportation[END_REF] We note that even if this result seems somewhat simple, it derives from non-trivial considerations such as the comparison of the critical distance and the area size. The important fact to retain here is that cost considerations in the case where a mode choice is available usually lead for large urban areas to leave the car and take the MRT.

Compiling data from 25 megacities in the world (see Appendix 1) for which we found an estimate of the population having access to the MRT (and therefore the quantity p), we compute the critical levels T * and P * (we note that for all these cities the definition of p is the same: it is the share of individuals living within 1km from a MRT station). For reasonable values of time [START_REF] Numbeo | [END_REF][START_REF] Crozet | Le temps et les transports de voyageurs (CEMT-ECMT)[END_REF], most of cities have T * = 0 and all have T * P (T * ranges from 0 to 50% of P ). Cases where T * = 0 mean that even at zero traffic the critical distance d(V, T = 0) is larger than the typical size of the city L. This depends a bit on the value of time but in most cases we do observe small values of T * indicating that cities are mostly in the saturated regime where the MRT is always more advantagous than the car. Public transports are so economical (compared to cars) that people living near rapid transit stations are highly likely to ride them. Thus, traffic does not appear as a determinant parameter in individual mobility choices as it concerns mostly individuals who have no choice but to drive and who suffer from onerous commuting costs and unavoidable time-consuming trips as traffic increases. We note here that correlations between the neighborhoods with low MRT density (small p) where individuals have expensive commuting costs and their revenue is an interesting field of study and appears as an important source of urban segregation [START_REF] Glaeser | Why do the poor live in cities? The role of public transportation[END_REF]. Since we have in general T * , P * P , we obtain the simple prediction that T P 1-p, a non-trivial consequence of rapid transit cheapness and individual choices of mobility. We compare the empirical car modal shares T P of these cities (see data description in Appendix 1) to our prediction on Figure 2 and observe -considering the simplicity of the model -a very good agreement and a relevant linear trend highlighting the efficiency of public transportation in reducing traffic. In particular, most of the European cities are well described by our prediction and we observe a few deviations. These discrepancies can probably find their origin in the existence of other modes of commuting, lower car ownership rates (e.g. in Buenos Aires [START_REF] Roque | An analysis of car ownership in Latin American cities: a perspective for future research[END_REF]), lower road capacities, higher cost of MRT, or a high degree of polycentrism. Discrepancies to the predicted value are probably mostly due to the existence of other modes of transport (walking or cycling), lower car ownership rates, or a higher cost of the MRT, etc.

Emitted CO2

Our model also provides a prediction for the transportrelated gas emissions and we will focus on the CO 2 case for which we obtained data. We make the simplest assumption where these emissions are proportional to the total time spent on roads. The quantity of CO 2 emitted for a driver residing at x is then given by Q

CO2 (x) ∝ d(x)[1 + (T /c) µ ] which leads to a total Q CO2 ∝ drivers i d(x i ) 1 + T c µ (11) ∝ g √ A(1 -p)P 1 + T c µ ( 12
)
where we used our result for the total traffic T = (1-p)P . We assumed that the sum i d(x i ) is of the form g √ AT where √ A sets the scale of displacement and where the prefactor g encodes the geometrical aspect of car mobility in the city, including the spatial distribution of residences and activities, and the transport infrastructure. Its estimation probably requires a more detailed, specific calculation but the important aspect here is the scaling with √ A (see also the appendices). We thus obtain that the annual CO 2 emitted by car and per capita is given by where τ = T c µ is the average delay due to congestion and is empirically accessible from the TomTom database [START_REF]TomTom Traffic Index[END_REF] (while the capacity c is usually not empirically accessible). It is interesting to note that Q CO2 is then the product of three main terms: size of the city × fraction of car drivers × congestion effects, which correspond indeed to the intuitive expectation about the main ingredients governing car traffic. Also, we note here that the dominant term of Q CO2 is proportional to the population indicating a simple linear scaling, and that nonlinear term could possibly appear as corrections, in contrast with previous results [START_REF] Louf | How congestion shapes cities: From mobility patterns to scaling[END_REF][START_REF] Bettencourt | Growth, innovation, scaling, and the pace of life in cities[END_REF], and which could explain the difficulties to reach a consensus about the behavior of this type of quantities (see for example [START_REF] Leitao | Is this scaling nonlinear[END_REF]). We compare our prediction Eq. 13 to disaggregated values of urban CO 2 emissions on Figure 3 and observe a good agreement. We observe some outliers like Buenos Aires which has a very small car ownership rate and thus lower than expected CO 2 emissions and areas such as New York which appears to be one of the largest transport CO 2 emitter in the world [START_REF]The OECD Metropolitan Areas Database visualized through the Metropolitan eXplorer[END_REF]. This result illustrates the role played by public transport and traffic in modulating transportrelated CO 2 emissions. Most importantly, we identify urban sprawl ( √ A) as a major criterion for transport emissions. We note that if we introduce the average population density ρ = P A , we can rewrite our result as

Q CO2 P ∝ √ A(1 -p)(1 + τ ) (13)
Q CO2 P ∝ √ P √ ρ (1 -p)(1 + τ ) (14)
i.e.

Q CO 2 P ∝ ρ -1 2
since √ P is a slowly-varying function within the scope of large urban areas. We understand here how Newman and Kenworthy [START_REF] Newman | Cities and automobile dependence: An international sourcebook[END_REF] could have ob-tained their result by assuming the density to be the control parameter. However, even if fitting data with a function of ρ is possible, our analysis shows that it is qualitatively wrong: the area size A and the public transport density p seem to be the true parameters controlling carrelated quantities such as CO 2 emissions. Mitigating the traffic is therefore not obtained by increasing the density but by reducing the area size and improving the public transport density. Increasing the population at fixed area would increase the emission of CO 2 (due to an increase of traffic congestion leading to an increase of τ ) in contrast with the naive Newman-Kenworthy assumption where increasing the density leads to a decrease of CO 2 emissions.

We also note that

∂ log Q CO 2 ∂ log p = -p 1
-p can be relatively large (in absolute value) while

∂ log Q CO 2 ∂ log A = 1
2 is small which suggests that the increase of transport density is in general more efficient than an area decrease (which is also less feasible). We also understand that if mass rapid transit is efficient in reducing transport-related greenhouse emissions, on the other hand reducing the road capacity c (to deter individuals from driving) may not be a good remedy, as we identified drivers with individuals with no other solution. Worse, by making their trips longer, diminishing c indirectly contributes to higher emissions and pollution as well as potentially socially segregating situations.

Average commuting time

Finally, we can also estime the average commuting time (details are given in Appendix 2) and obtain for the one-trip commuting time τ c averaged over the population following expression

τ c = p f + g √ A v m + (1 -p) g √ A v c (1 + τ ) ( 15
)
where g is a geographical factor that encodes the spatial complexity of trips. The comparison with the data displays relatively large fluctuations, but our analysis seems to capture the main trend and the one-parameter fit over g using our expression Eq. 16 leads to the average value g ≈ 0.203. We note that for a uniform distribution of residences, a simple calculation leads to g = 2/3 √ π ≈ 0.376 and in the simple isotropic case where the density decreases with the distance r to the center as ρ(r) = ρ 0 (1r/L), we obtain g = 1/3 √ π ≈ 0.188. The average value g obtained by the fit is in between these two theoretical estimates. However, polycentrism and more generally the spatial organization of the city and its infrastructure certainly play an important role here and our model can only provide a first approximation to the average commuting time.

DISCUSSION

We presented a simple model for the car traffic and its consequences in large cities. This model illustrates how a combination of statistical physics, economical ingredients and empirical validation can lead to a robust understanding of systems as complex as cities. In particular, this approach is in contrast with previous works where it was a priori assumed that urban density was pivotal, and our aim here was to capture the essence of the urban mobility phenomenon. Our analysis shows that traffic related quantities are governed by three factors: access to mass rapid transit, congestion effects and the urban area size. In order to reduce CO 2 emissions for example, our model suggests to increase public transport access either by increasing the density around MRT stations or to increase the density of public transport (in contrast with the conclusions of an econometric study in the US [START_REF] Duranton | The fundamental law of road congestion: Evidence from US cities[END_REF]), or to reduce the urban area size (impossible in most contexts). Increasing the cost of car use seems actually unable to lower car traffic in the absence of alternative transportation means. Also, it is important to note that increasing in general the density at fixed both area size and MRT access would actually increase CO 2 emissions, in sharp contrast with the commonly accepted paradigm about the effect of urban density. Finally, we insist on the fact that in this model we voluntarily left out a number of parameters such as other transportation modes (buses, tramway), polycentrism, the transport network structure, fuel price and tax, dynamic road pricing, etc. but our main point was to fill a gap for understanding traffic in urban areas by proposing a parsimonious model with the smallest number of parameters and the largest number of predictions in agreement with data. Given the simplicity of this model we cannot expect a perfect agreement with data for various and different cities, but it seems that this approach captures correctly all the trends and identifies correctly the critical factors for car traffic. This seems to be a basic requirement before adding other factors and increasing the complexity of the model. Also, it seems at this stage necessary to first encourage the measure and sharing of data such as the density of public transport in order to propose further tests of our theoretical framework.

• The definition of a metropolitan area varies from one country to another but we aimed at assembling a statistically coherent set of agglomerations. Populations and areas were collated from Wikipedia data on metropolitan statistical areas in accordance with national definitions.

• Various indicators were compiled from diverse sources: the quantity p is taken from transit accessibility reports [START_REF] Marks | People Near Transit: Improving Accessibility and Rapid Transit Coverage in Large Cities[END_REF][START_REF] Singer | Fast Cities. A comparison of rapid transit in major Canadian cities[END_REF][START_REF]Part de la population à proximité d'un arret de métro ou d'un tram chrono[END_REF][START_REF]Singapore's Transportation Challenges and Plans[END_REF][37]41]. A common metric is the share of individuals living within 1 km from a MRT station (this is the assessed maximal walking distance). We note however that this is an unusual indicator of mobility that is not always easy to find and was here the main bottleneck for getting a larger set of cities. Modal shares were obtained from [START_REF] Mckenzie | Who drives to work? Commuting by automobile in the United States: 2013[END_REF][START_REF]Commuters using sustainable transportation in census metropolitan areas Census of Population 2016[END_REF][START_REF][END_REF][41][START_REF]EMTA Barometer of Public Transport in European Metropolitean Areas[END_REF][START_REF] Atm | Ajuntament de Barcelona, Enquesta de Mobilitat en Dia Feiner[END_REF][44], CO 2 emissions from [START_REF]The OECD Metropolitan Areas Database visualized through the Metropolitan eXplorer[END_REF][46][START_REF]Cities and climate change: an urgent agenda[END_REF][START_REF]GHC Inventory mitigation[END_REF][49], and commuting time from [START_REF] Lebrun | Cahiers de lObservatoire de la mobilité de la Rgion de Bruxelles-Capitale, Les pratiques de déplacement à Bruxelles[END_REF][START_REF]Census of Population[END_REF][START_REF] Uranga | Los portenos pierden hasta 18 dias al ano en ir y volver de trabajar[END_REF][START_REF] Pereira | differences between metropolitan areas, by income levels and gender[END_REF][START_REF]Governing calculations of 2011-2015 American Community Survey data[END_REF][START_REF] Bitre | Five facts about commuting in Australia[END_REF][56][START_REF]Statistics Korea, Sample Enumeration Results of the 2005 Population Census (Internal migration and commuting)[END_REF][START_REF]Beijing workers have longest daily commute in China at 52[END_REF].

• Values of time were assessed by taking half of the hourly wage after tax for each city [START_REF] Numbeo | [END_REF][START_REF] Crozet | Le temps et les transports de voyageurs (CEMT-ECMT)[END_REF].

• Congestion delays were taken from the TomTom index [START_REF]TomTom Traffic Index[END_REF], except for Seoul and Tokyo where we had no data (we assumed an average delay of 50%).

• The road capacity was computed from the congestion delay τ by c ∼ P/τ 1/µ with the value µ = 2 that we used for all cities.

• For the velocities v m and v c and the costs C c and f we choose the same values for all cities. For v m wikipedia data [62] displays values in the range 25 -35 km/h depending on the city and we took v m ≈ 30km/h. The free flow car velocity v c also depends a bit on the city and varies from 30 km/h (without congestion effect) in European cities such as Paris to 56 km/h in some american cities [63]. We decided to take an average value v c ≈ 40 km/h.

• For C c , we used a cost simulator [64] which gives on average a value of C c ≈ 15 USD per trip.

• For f we counted on average 10 minutes to reach the station, 10 minutes to reach the office and 10 minutes for transit for a total of about f ≈ 30 minutes. This value can certainly be improved but it is difficult to get and we expect it not to vary too much (at least not over more than one order of magnitude).

Appendix 2: Average commuting time

The commuting time for MRT users is given by f + g √ A/v c while for car users it is g √ A(1 + τ )/v c where g is a geographical factor. We therefore obtain for the onetrip commuting time τ c averaged over the population the [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF]. We perform the fit on the single parameter g using Eq. 16, leading to an average value g ≈ 0.203 (the Pearson correlation is here 0.65).

following expression

τ c = p f + g √ A v m + (1 -p) g √ A v c (1 + τ ) (16)
The comparison of this result with empirical data is shown on Figure 4. Even if we observe relatively large fluctuations, our analysis seems to capture the main trend and the one-parameter fit using our expression Eq. 16 leads to the average value g ≈ 0.203. We note that for a uniform distribution of residences, a simple calculation leads to g = 2/3 √ π ≈ 0.376 and in the simple isotropic case where the density decreases with the distance r to the center as ρ(r) = ρ 0 (1r/L), we obtain g = 1/3 √ π ≈ 0.188. The average value g obtained by the fit is in between these two theoretical estimates. From Fig. 4 we can compute the 'effective g eff ' for each city and compare it to the average value g. This quantity g eff (and the ratio η = g eff /g) encodes both the complexity of the population and activity densities, and of the transportation infrastructure. For our set of cities, the average ratio is of order 1.37 and we observe outliers with small values such as Barcelona (η = 0.04), Seoul (η = 0.62), and large ones (Buenos Aires 1.92, Rotterdam 2.32, Singapore 3.38). Most of the cities (76%) have however here a ratio η larger than 1 which implies that our model underestimates in general the commuting distance. Many local effects can explain the variations observed in the commuting distance: the existence of polycentricity could reduce the commuting distance (for example, it seems that we overestimate the commuting time for New York which might be a consequence of polycentrism), but other factors such as poor transportation infrastructures (or traffic bottlenecks due to geographical constraints) could have the opposite effect of increasing this quantity. In general we expect that the transport network structure will have an important impact, especially if it is very anisotropic in space. The determination of the commuting distance for each city probably implies to take into account a large number of specific details, but our analysis shows that it is consistent to assume that the order of magnitude of this quantity is √ A.

FIG. 1 :

 1 FIG.1:The most advantageous mode of transport depends on the value of time of individuals and the distance to the urban center. The limit between the two areas evolves with congestion: the larger the traffic (curves from blue to green) and the larger the area in which rapid transit is beneficial compared to car driving. The grey solid vertical line corresponds to the size of the urban area and indicates the critical value of time (dashed red line) below which rapid transit is advantageous in the whole agglomeration whatever the value of congestion. The values of the parameters are chosen here as: Cc = 15 $, vc = 40 km/h, vm = 30 km/h, f = 30 minutes (see the Appendix 1 for a precise description of all data).

FIG. 2 :

 2 FIG.2: Comparison between the observed car modal share T /P and the share of population p living near rapid transit stations (less than 1 km) for 25 metropolitan areas in the world. The red line is the prediction of our model (R 2 = 0.69). Discrepancies to the predicted value are probably mostly due to the existence of other modes of transport (walking or cycling), lower car ownership rates, or a higher cost of the MRT, etc.

FIG. 3 :

 3 FIG.3: Comparison between the annual transport-related CO2 emissions per capita and the effects of congestion, area size and rapid transit density predicted by our model. The red line is the linear fit of the predicted form y = ax where a ≈ 0.064 CO2tons/km/hab/year (the Pearson coefficient is 0.79) We had no congestion estimate for Seoul and Tokyo and we used an average congestion rate τ = 50%.

FIG. 4 :

 4 FIG.4: Average commuting time τc measured for different cities versus the predicted value of our model (Eq. 16). We perform the fit on the single parameter g using Eq. 16, leading to an average value g ≈ 0.203 (the Pearson correlation is here 0.65).

Acknowledgments. VV thanks the IPhT for support during his internship.

Appendix 1: Data sources

We studied 25 metropolitan areas from Europe, America, Asia and Australia. The number of cities was limited by the availability of data on MRT accessibility and reliable modal share estimates. All the data is freely available and we list here their sources.


Andreu Arinyo-I-Prats

Lluís Mas-Ribas

Jordi Miralda-Escudé

Ignasi Pérez-Ràfols

Pasquier Noterdaeme

A metal-line strength indicator for Damped Lyman Alpha (DLA) systems at low signal-to-noise

, showing that DLAs are associated with relatively massive halos. However, the low resolution and signal-to-noise of BOSS spectra do not allow precise measurements of the DLA metal lines. We define a metal strength parameter, S, based on combining equivalent widths of 17 metal lines, that can be measured with an optimal signal-to-noise ratio for individual DLAs in BOSS spectra, allowing for the classification of the DLA population into subgroups of different S. We present the distribution of this DLA metal strength and the dependence of its mean value on N HI and redshift. We search for systematic effects and variations in the catalogue purity by examining the dependence of the S distribution on the spectral signal-tonoise and the estimated error on S. A catalogue of DLAs with measured equivalent widths for the selected 17 metal lines and the value of S is made publicly available, which will be used to measure the dependence of the DLA bias factor on the S parameter. The relation of the metal strength on the gas metal abundances and velocity dispersion can be constrained by studying the stacked metal absorption spectra of DLAs as a function of S, allowing for future determinations of the dependence of the bias factor on the metallicity and velocity dispersion of DLAs.

INTRODUCTION

Observations of Damped Lyman Alpha systems (DLAs) in quasar absorption spectra are crucial to understand galaxy formation at different epochs (e.g., Wolfe et al. (1986Wolfe et al. ( , 2005))). Most of the atomic hydrogen in the Universe is present in DLAs, or absorption systems with N HI > 2 × 10 20 cm -2 , which are self-shielded against the external ionizing background. This atomic hydrogen accounts for ∼ 2 to 3% of all the baryonic matter in the Universe, a fraction that remains roughly constant during the epoch of the maximum rate of star formation in the Universe at 2 < z < 4 Noterdaeme et al. (2012) and is comparable to the fraction of baryons in stars at that epoch Shapley (2011).

A moderately high large-scale bias factor for DLAs of b DLA 2 has been measured by Font-Ribera et al. (2012); Pérez-Ràfols et al. (2018) at z 2.3, implying that most DLAs are hosted by halos of massive galaxies or galaxy groups, and that an extended distribution of atomic gas is generally present in these halos out to radii larger than the size of visible galaxies to account for the large cross sections required to explain the observed rate of incidence. Models of galaxy formation require strong winds to explain this high bias factor, which can expel the gas from low mass halos to reduce the contribution of dwarf galaxies to DLAs, and spread the gas out to large radius in massive halos to increase their DLA cross section (e.g., Barnes & Haehnelt (2014); Bird et al. (2015)).

Although the Lyα absorption profile of a DLA tells us only the hydrogen column density, much more information can be obtained from the DLA metal lines: the abundances are typically similar to those of metal-poor globular clusters and halo stars in the Milky Way, and the absorption line profiles reveal a rich and diverse velocity structure characterized by multiple components, with a total velocity width ranging from 10 km s -1 to 200 km s -1 (Prochaska & Wolfe 1997, 1998;Wolfe 1998). The derived metallicities are distributed over a broad range (10 -3 ≤ Z/Z ≤ 10 -1 ), and their average declines slowly with redshift (e.g., Prochaska & Wolfe (2002); Rafelski et al. (2012aRafelski et al. ( , 2014))).

Most studies of the DLA metal lines are based on highresolution, high signal-to-noise quasar spectra, which are necessary to detect the small equivalent width lines and to reveal the smallscale complexity of the velocity profiles. However, the large samples of DLAs required to measure cross-correlations that are used to infer the DLA linear bias factor can only be obtained at present with low resolution and low signal-to-noise spectra similar to those of the Baryon Oscillation Spectroscopic Survey (BOSS) of SDSS-III Eisenstein et al. (2011); Dawson et al. (2013); Smee et al. (2013); Alam et al. (2015).

Our aim in this work is to define a parameter characterizing the strength of the metal lines of a DLA, which can be measured with the highest possible signal-to-noise ratio by combining all the metal absorption lines that are usually observable for DLAs. If this parameter can be measured sufficiently accurately even in low signal-to-noise spectra similar to BOSS, then the DLAs can be classified into groups of different metal strength, for which several average properties such as the large-scale bias factor can be measured as a function of the metal strength. We may hope in this way to test, for example, the presence of a relation between host halo mass and gas metallicity in a DLA, which is naturally expected because of the well-known mass-metallicity relation for galaxies. As we shall explain below, caution is required to reach any such interpretations because the metal strength parameter we will define depends only on the metal line equivalent widths that are measurable in low-resolution absorption spectra, which are affected both by metal abundances and the gas velocity dispersion. Nevertheless, the mean absorption spectrum can also be measured for groups of DLAs with different metal strength using the technique discussed in Mas-Ribas et al. (2017), which contains information that can help disentangle variations of the gas velocity dispersion and metal abundances with metal strength.

This work is organized as follows: §2 specifies the DLA catalogue and the quasar spectra from BOSS that we use. In §3 we describe in detail our method to measure line equivalent widths for a set of 17 low-ionization metal lines that have the highest equivalent widths in DLAs. In §4 we define the metal strength, and a second quantity that is corrected for the effect of the hydrogen column density on the mean strength of metal lines. Results are presented in §5, where we discuss the distribution of this metal strength parameter and several systematic effects due to impurities in the DLA catalogue, and we make publicly available a catalogue with our measurements of equivalent widths and the metal strength parameter for each DLA. Finally, we present the conclusions in §6. Applications to measure the mean bias factor and other properties as a function of metal strength will be presented in future papers.

DATA SAMPLE

We use the quasar spectra in the complete SDSS-III BOSS Data Release 12 (DR12), from the Quasar Catalogue DR12Q Pâris et al. (2017). A detailed description of the SDSS telescope and the BOSS instrument obtaining the spectra is found in Gunn et al. (1998Gunn et al. ( , 2006)); Dawson et al. (2013); Smee et al. (2013), and the method to select quasar targets is described in Ross et al. (2012).

We use the DR12 extension of the DLA catalogue of Noterdaeme et al. (2012), containing a total of 34, 050 DLA candidates with a column density log(N HI ) ≥ 20. The method to detect these DLAs, described in Noterdaeme et al. (2009), is an automatic profile recognition procedure using Spearman correlation analysis with a Voigt profile. Only the Lyα absorption line is considered to decide if an absorption feature is included as a DLA in the catalogue. This ensures that there is no selection bias in favor of DLAs with strong metal lines. The presence of metal lines is, however, used to refine the accuracy of the measured redshift of the detected DLAs. This catalogue will be referred to as DR12-DLA from now on.

The DLA sample we use in this work is very similar to the sample of Mas-Ribas et al. (2017) (who used an earlier version of the DR12-DLA catalogue with minor differences in the selected DLAs), who showed the distribution in redshift and N HI in their figure 1. Most DLAs are at redshifts 1.9 < z < 3.5, with a small fraction extending to higher redshifts. Although DLAs are usually defined as having log(N HI ) > 20.3 (roughly the column density at which the gas becomes mostly neutral due to self-shielding, although this depends on the gas density; see Wolfe et al. (1986)), we use all the systems going down to column densities N HI > 10 20 cm -2 to have a larger sample.

We note that these systems are only DLA candidates, and that some fraction of them are expected to be false DLAs arising from a combination of noise and regions of strong Lyα forest absorption that are confused with a DLA in low signal-to-noise spectra. The purity of the catalogue, or fraction of DLAs that are real, is expected to decline as the signal-to-noise ratio and the column density decrease.

The catalogue gives, for each detected DLA, the quasar and DLA redshift, the continuum-to-noise ratio (hereafter, CNR) of the spectrum in the Lyα forest region (as defined in Noterdaeme et al. (2009)), and the DLA column density. For each DLA in the catalogue, we use the corresponding quasar spectrum from BOSS to measure the associated metal lines. We use the co-added spectra, which are provided in wavelength bins of width ∆ log 10 (λ) = 10 -4 , corresponding to a velocity width of 69.05 km s -1 .

MEASURING METAL LINE EQUIVALENT WIDTHS

In this section we describe the method used to compute equivalent widths (W) of a selection of 17 metal lines of DLAs. A weighted average of these equivalent widths will be used as a definition of the metal strength of each DLA. Briefly, our method uses a measurement window around the central wavelength of the absorption line over which the equivalent width is integrated, and two windows around it to determine a continuum from a linear regression of the measured flux. In §3.1 we explain the selection of the 17 lines. The windows are described in §3.2. In §3.3 we specify which of the 17 lines are used for each individual DLA depending on how the windows are placed in the quasar spectrum, and §3.4 gives the details of how the equivalent widths and their errors are calculated.

The DR12-DLA catalogue of Noterdaeme et al. (2009) includes measurements of W for 10 metal lines. However, their method provides a biased estimate of these equivalent widths because they are only measured when their detection is considered significant. This means that a metal line may be included in the catalogue when the value of W has been increased by noise, or may be dropped when it has been reduced, systematically affecting the average W. Negative equivalent widths caused by noise also need to be included in the catalogue to avoid bias. Our method will select the lines to be measured depending only on the location of the windows to be used for the continuum determination and the equivalent width integration, but not on the value of W that is derived.

List of Metal lines included

The 17 metal lines we use to evaluate our metal strength for each DLA are selected to be transitions of low-ionization species that lie in the DLA rest-frame wavelength interval between 1260 Å and 3000 Å, with a mean equivalent width measured from the stacked Prochaska et al. 2001). The first three lines are blends of two lines and we give the wavelength of the strongest line, except for O I-Si II (which are comparably strong) where we give the average of the two wavelengths. In the third blend, C II * is a metastable state of CII. Third column: Equivalent width measured from the stacked spectrum in (Mas-Ribas et al. 2017). Fourth column: Mean equivalent width W calculated with our method as described in §3 for spectra with CNR > 2 as in (Mas-Ribas et al. 2017). Fifth column: Mean contribution C k of each line to the total metal strength in the DLAs where W of the line is measured (see §4 for details). Sixth column: Fraction x k of the 34, 050 DLAs in which each line is actually measured. Seventh column: Number of DLAs in which each line has been measured. The low-ionization species are either neutral atoms or ions that are once or twice ionized, but we exclude higher ionization lines like C IV and Si IV. We use this combination of low-ionization transitions to obtain a metal strength that reflects a property of the low-ionization gas with the highest possible signal-to-noise, by combining W measurements of all the available lines. We exclude high-ionization lines because these are known to reflect a more extended gas distribution with different physical properties. In practice, we know that this metal strength will depend both on the metal abundances and velocity dispersion, because of line saturation effects.

Name

λ [ Å] Stacked W [ Å] W [ Å] C k x k N k Si II-
The 17 selected line transitions are listed in Table 1, where we give the equivalent width computed from the stacked DR12 spectrum designated as total in Mas-Ribas et al. ( 2017), compared to the mean equivalent width W measured with our method, as described below. The first three lines in Table 1 are blends of lines that are too close to be measured separately, and we list the total equivalent widths for the blend. Their wavelengths are that of the Si II transition in the first blend (which is much stronger than Fe II), the average of the O I and Si II transitions in the second blend (which are comparably strong), and the C II transition in the third blend (which is much stronger than the transition of the metastable state C II * ). The rest of the quantities will be discussed in §4.

Measurement and Continuum Windows

We use a fixed measurement window to integrate the equivalent width for all the lines, with a width that is normally set to 15 pixels in the co-added spectra, except for the O I-Si II blend for which we use a window width of 27 pixels. The measurement window is centered on the pixel that includes the central wavelength of the line listed in Table 1 multiplied by 1 + z d , where z d is the DLA redshift, and includes 7 pixels on each side (13 for O I-Si II). The width of 15 pixels, corresponding to a velocity width 1036 km s -1 , is generally wide enough to include most of the absorption components in DLAs, as can be seen from the metal line profiles in the stacked absorption spectrum of DLAs in Mas-Ribas et al. (2017) (see their Figure 10), which have a dispersion of ∼ 100 to 150 km s -1 (only moderately wider than the Point Spread Function of the BOSS spectrograph wavelength resolution). However, the equivalent widths we measure are underestimated if there are large DLA redshift errors in the DLA catalogue we use, which shift the metal lines partly outside the measurement window. In addition, the O I-Si II blended line requires the wider window mentioned above because the two lines are resolved and are spread over a wider interval than the other blends. We will show below the impact on the mean equivalent widths of decreasing the measurement window width to 10 pixels or increasing it to 20 or 25 pixels (Figure 2 and§5). For these cases, we also change the O I-Si II blend to 25, 29 and 31 pixels, respectively.

The windows for determining the continuum are both set to a width of 21 pixels, on the left and the right of the measurement window. We generally leave 1 pixel that is not used between the end of the measurement window and the start of each continuum window, with the exception of some lines for which the continuum window placed in this way would include another DLA metal line that would systematically lower the estimated continuum. These exceptions are listed in Table 2, where we give the space left between the center of the line and the start of the two continuum windows on each side for these set of lines (for all other lines, the normal space left is 8.5 pixels, corresponding to half of the 15 pixel width of the measurement window plus one). The size of the two continuum windows remains fixed at 21 pixels for all these lines.

Table 2. Range around each metal absorption line at which the two continuum windows start, given in pixel numbers (p) and velocity interval (s). The two continuum windows of any of the two lines in each of the 4 pairs start a number p of pixels to the right and to the left of the line center. When measuring one line in any of these four pairs, the other line is avoided by leaving it in between the integrating and the continuum windows.

Excluded metal lines for each DLA

DLA metal lines often fall in regions of the spectrum where the equivalent width measurement is subject to large systematic errors, due to the presence of quasar emission lines or of sky lines from the atmosphere that increase the noise. To reduce these systematics we measure only a subset of the 17 metal lines for each DLA that satisfy the following criteria:

• The right end of the continuum window (at the longest wavelength) must be below the longest wavelength of the BOSS spectrum.

• The left end of the continuum window must lie at least 30 pixels redwards of the center of the Lyα emission line of the quasar (at λ = 1215.6 Å), to avoid the Lyα forest and the effect of the Lyα quasar emission line on the measurement of W.

• The left and right ends of the continuum windows must be separated by more than 30 pixels (or 2071.5 km s -1 ) from the center of any of the following four quasar emission lines: Si IV-O IV at λ = 1400.0 Å, C IV at λ = 1549.2 Å, C III at λ = 1908.7 Å, and Mg II at λ = 2798.7 Å. These are the strongest quasar emission lines (see, e.g., Table 3 in Pâris et al. (2012)), which we avoid because they affect the continuum determination and the measurement of W.

• The measurement window must not include any skyline. We discard the metal line whenever the central wavelength of any skyline is within the measurement window or in one of the two adjacent pixels, to avoid the increased systematic error on the equivalent width caused by skylines. Following previous analyses of BOSS quasar spectra, we use the set of 872 skylines 1 .

Calculation of the metal line equivalent width and error

The quasar continuum around a metal absorption is estimated by doing a linear regression of the flux values f i in every pixel i in the two continuum windows to the left and right of the measurement window, which contain a total of 42 pixels. First a linear regression 2 is computed applying the inverse variance w i = 1/σ 2 i of the pixels i as weights, where σ i is the noise. This yields a preliminary determination of a continuum c i , from which we compute the transmission fraction F pi = f i /C pi . Next, we eliminate outliers in the continuum windows, which may result from highly noisy flux measurements in skylines. Note that DLA metal lines are discarded only when 1 https://github.com/igmhub/picca/blob/master/etc/ dr14-line-sky-mask.txt 2 http://www.statsmodels.org/stable/generated/ statsmodels.regression.linear_model.WLS.html skylines are present in the measurement window, but not in the continuum window. We eliminate any pixels for which |F pi -1| > 4 as outliers. We repeat the determination of the continuum and we check again for outliers, until none are left. If this process leaves less than 5 pixels in any of the two continuum windows, the line is discarded and not measured for the DLA being analysed.

Then we repeat the determination of the continuum with nonweighted linear regression with the pixels that are left after elimination of outliers. The reason for not using weights in our final determination of the continuum is that the weights that minimize the impact of more noisy pixels also introduce a bias that systematically lowers the continuum, because pixels with less flux are assigned lower noise owing to the photon noise contribution calculated by the BOSS pipeline. This systematic effect causes equivalent widths to be underestimated due to the lowering of the continuum, by an amount that increases with spectral noise and introduces an unwanted artificial correlation of our metal strength parameter and its error, so we decided not to use weights to obtain the final continuum. With the final continuum estimation C i , we compute the transmission fraction F i = f i /C i . At the same time, not using the weights for our final continuum implies that outliers caused by skylines or other large errors can strongly distort the determination of the continuum and they need to be eliminated.

We checked that different ways of computing the linear regression, weighted or unweighted, and varying the outlier condition to |F pi -1| > (2, 3, 5), have only percent effects on the mean equivalent widths.

Figure 1 Finally, we obtain the rest-frame equivalent width of a line k simply by summing the absorbed fraction over all the N m pixels of the measurement window,

W k = N m i=1 (1 -F i )λ k p 0 , (1)
where λ k is the central rest-frame wavelength of the metal line being measured, and p 0 = 10 -4 × log(10) = 2.303 × 10 -4 is the width of the BOSS pixel in log λ. Note that some pixels in the measurement window may have F i > 1, and in fact, some equivalent widths we compute are negative because of the noise. Still, as explained before, they need to be included to have unbiased mean properties. The error on this equivalent width is computed as

k = σ F N 1/2 m λp 0 , (2)
where σ 2 F is the average of the squared transmission error, σ 2 Fi = σ 2 i /C 2 i , of the pixels in the two continuum windows that are used to determine the continuum linear regression. We use the noise in the continuum windows (instead of the measurement window) to estimate the error k because of the correlation of σ Fi and F i discussed earlier, which can induce an artificial correlation of the estimated value of W and its error. We note that the error k accounts only for pixel noise, assuming that it is uncorrelated in all the pixels, but excludes any systematic errors due to the continuum fitting.

We have tested our method by inserting mock lines in empty parts of a spectrum and successfully recovering their equivalent widths within the estimated uncertainty. We checked that equivalent widths estimated in spectral intervals where we do not expect any absorbing line return measurements consistent with zero.

As a final test, in Figure 2 we compare the mean equivalent width we derive for each of the 17 lines in Table 1 with the value obtained from the stacked absorption spectrum presented in Mas-Ribas et al. ( 2017), and we also check the sensitivity of this mean equivalent width to the size of our measurement window by varying it from our standard width of 15 pixels to 10 and to 20 pixels. The mean values generally agree. The largest discrepancy occurs for the Fe II 2344 line, for which our method yields an equivalent width ∼ 20% smaller than in Mas-Ribas et al. (2017). Our derived W generally increases with the measurement window width, as expected because of the effect of redshift errors. The overall agreement is relatively good for a 15 pixel width, which we choose as the optimal and standard one from this plot.

DEFINITION OF THE METAL STRENGTH

We now define a metal strength for each DLA intended to provide a weighted average of the strength of the observable associated metal lines with an optimal signal-to-noise, which can be used in spectra with high noise where individual lines are generally barely detectable. If the metal lines that have been measured for a given DLA are labeled by the index k, with equivalent widths W k and error k from equations ( 1) and ( 2), the metal strength is defined as

S = k W k / k 2 • W k /W k k W k / k 2 , ( 3
)
where W k is the mean equivalent width of line k over all the DLAs for which it is measured. The metal strength is equal to unity if all the lines of a DLA have the mean value W k , and is in general a measure of the ratio of the equivalent widths to their mean value. Each line is weighted by (W k / k ) 2 , which is the expected squared signal-to-noise ratio of the measurement if the line has an equivalent width equal to the mean. These mean values are listed in the fourth column of Table 1. The error of this metal strength is computed from the errors of each individual line as

S =       k W k k 2      -1/2 . ( 4
)
This metal strength parameter reflects how strong the metal lines of a DLA are, and is chosen as a quantity that can be measured with an optimal signal-to-noise ratio from spectra similar to those in the BOSS survey. The metal strength is expected to increase both with the metal abundances and the velocity dispersion. A larger velocity dispersion implies a wider line which is less saturated, thereby yielding a larger equivalent width for a fixed column density of metals.

To visualize the metal lines that are most important to determine S, we define the mean contribution C k of metal line k as follows:

C k = 1 N k N k i=0 W k / k,i 2 j W j / j,i 2 , (5)
where the sum over i includes all the N k DLAs in which the metal line k has been measured, j,i is the error of the metal line j in the DLA i, and the sum in the denominator adds the weights of all metal lines j that have been measured in the DLA i. This means that C k is the mean contribution to the metal strength S from line k, restricted only to the fraction x k of DLAs in which this metal line has been measured. The average contribution to all the DLA sample is the product C k x k , which adds to unity when summed over the 17 metal lines. These quantities are listed in Table 1, and shown also in Figure 3. The lines that are important for most DLAs, with the highest value of C k x k , are the first 6 at short wavelength. A few of the longer wavelength ones are important only in a small fraction of DLAs. The fraction x k is small for the first line, the SiII-FeII blend, which has the shortest wavelength, because any lines that are close to the quasar Lyα emission line or that fall in the Lyα forest region are excluded. Then, x k increases with increasing line wavelength until it starts decreasing rapidly above ∼ 2000 Å due to the exclusion of metal lines when a sky line falls within the measurement window, and the large abundance of sky lines in the red part of the spectrum.

The distribution of the number of lines contributing to the

Correction for the dependence of W k on N HI .

We expect the mean equivalent width of any metal line to increase with N HI , because for a fixed metal abundance and internal velocity distribution, the metal column density should be proportional to the hydrogen column density. This dependence was confirmed in Mas-Ribas et al. ( 2017) from stacked spectra of DLAs in different bins of N HI (see their Table 3). Our goal is to classify DLAs depending on a parameter that reflects only their metal content and velocity For this purpose, we define the HI-corrected metal strength, S HI , in the same way as S in equation 3, but replacing the constant values of W k by mean equivalent widths that depend linearly on log(N HI ):

W k (N HI ) = a k (log(N HI ) -20) + b k . ( 6
)
We fit the parameters a k and b k by dividing the DLA sample into the following 5 intervals of log(N HI ): [20.0, 20.20), [20.2, 20.4), [20.4, 20.63), [20.63, 21.0), and [21.0, 22.2]. Then, for each metal line, we compute the mean value of W k and log(N HI ) in the five intervals, and we make an unweighted linear regression of these 5 points.

The result is shown in Table 3, in columns 4 and 5. These linear regressions to the mean equivalent width as a function of N HI were also obtained in Mas-Ribas et al. (2017) from stacked spectra, and their results are shown in columns 2 and 3. The values are in general similar, and differences can be attributed to the different methods of determining the quasar continuum and weighting the DLAs.

The mean equivalent widths calculated from equation 6 for the column density of each DLA are then used to compute the corrected metal strength S HI with the same equation 3 as for S.

RESULTS

We now present the distribution of the metal strength parameter S and its variation with column density and redshift. We also examine if the metal strength distribution has any dependence on the continuum-to-noise ratio in the Lyα forest, and the computed error S , as a test of any possible variation of the purity of the DLA catalogue that might affect derived statistical properties of DLAs as a function of S. As the main product of this work, we present the catalogue of metal line equivalent widths and the metal strength parameter with its error, both corrected and uncorrected for its N HI dependence, for all the DLAs in DR12-DLA.

The metal strength distribution

The distribution of the metal strength S is shown in the top panel of Figure 5, including only DLAs with an error S < 0.5. In general, we exclude DLAs with larger error to better see their intrinsic properties. The distribution we measure is the convolution of the true distribution with the error distribution. Systems with negative S should be a consequence of the error, and this is consistent with the shape of the curves at S < 0, except at S < -1 where the extended tail indicates the presence of non-Gaussian errors that are likely due to systematics in the continuum fitting. The S distribution becomes wider with the measurement window width owing to increased errors. The N HI correction practically does not affect this distribution.

The bottom panel shows the distribution of S . The error increases with the measurement window simply because of the increased number of pixels in equation 2. The majority of DLAs have errors S < 0.5, and selecting this subset helps make the distribution of S we obtain in the top panel less distorted from the true one. We see that this subset of DLAs have an error S smaller than the true dispersion in S, and therefore can be divided into smaller subsets corresponding to intervals in S to measure other DLA properties as a function of S. In future papers, we will study the dependence of the DLA bias factor and the mean absorption spectrum on the parameter S.

Results in the rest of the paper are shown only for DLAs with S < 0.5, unless otherwise specified.

Metal strength dependence on N HI and redshift

Figure 6 shows the mean value of S in several column density intervals as the green points. As expected, S increases with N HI simply due to the increasing metal column densities. After applying the correction discussed in §4.1 (blue points), the mean value of S HI becomes indeed nearly constant. Note that the mean value of S should be unity for the whole DLA sample; when we eliminate DLAs with S > 0.5, the mean value of S increases slightly above unity because of a decrease of the mean S with S that will be discussed below. in Figure 7. The mean S is roughly constant at z < 3.4, and then drops rapidly at higher redshift. Some of this decline should be due to the well-known decrease of the average DLA metallicity with redshift (e.g., Kulkarni et al. (2005); Rafelski et al. (2012b)).

A decreasing velocity dispersion with redshift may also contribute because the mean value of S increases with velocity dispersion at a fixed metal column density owing to the effects of line saturation. However, the dependence of the purity of the catalogue on redshift may be the main effect causing the rapid drop at z > 3.4. Only ∼ 10% of our DLA sample is at z > 3.5 (see Fig. 1 of Mas-Ribas et al. ( 2017)), and the difficulty in detecting DLAs is greatly increased at high redshift because of the increased mean absorption of the Lyα forest, which is likely decreasing the purity of the catalogue. More complete studies revising the DLA detection method will be needed to resolve this issue.

Dependence of the S distribution on S and the purity of the DLA catalogue

Our main goal in defining and measuring the metal strength S for individual DLAs is to be able to measure average quantities of the DLAs as a function of S. An important difficulty to accomplish this is the imperfect purity of the catalogue. For example, false DLAs arising purely from spectral noise should have a mean value of S and bias factor equal to zero, and this can induce a spurious increase of the bias factor with S. Other false DLAs may be dense Lyα forest regions with a column density and mean S much lower than DLAs, and a different bias factor.

A first test we can do to investigate possible variations of the purity with S is to check if the mean value of S depends on the continuum-to-noise (CNR) ratio of the Lyα forest region, where the DLA absorbers are searched. We expect the purity of the catalogue to increase with CNR. This is shown in the top panel of Figure 8, both for the uncorrected and N HI -corrected cases. The mean S is consistent with a constant value for CN R > 3, suggesting that as long as we restrict any analysis to this minimum CNR, there should not be strong effects of the catalogue impurity caused by spectral noise. The bottom panel shows that, as expected, the error S decreases with CNR (since they both reflect the noise amplitude in the spectrum, although on different wavelength ranges). In both the uncorrected and corrected cases, the sample is restricted to DLAs with S < 0.5, or SHI < 0.5.

A second test we do is to examine if the S distribution varies with the error S . In principle, the true distribution of S should not depend on S , which is determined using only the flux measured in spectral pixels in the continuum windows and is independent of the flux values in the measurement window, as explained in §3.4. However, in reality there is a fairly strong dependence as we show in Figure 9, where we have divided DLAs into three intervals of S ([0.0, 0.16), [0.16, 0.33), and [0.33, 0.5]). If the true S distribution stayed constant, we would expect only an increased width of the distribution as S increases because of the convolution with the error, but we also observe a decrease of the mean value of S with increasing S . The effect is nearly the same for the corrected S HI , as seen in the bottom panel. This effect is more clearly seen by examining the mean value of S as a function of S , shown in Figure 10. Our standard upper limit S < 0.5 is shown as the dashed vertical bar, and we divide the range 0 < S < 0.5 into 5 intervals. Blue points in the top panel show that in our standard case, the mean value of S declines from 1.2 to 0.8 as S increases up to its upper threshold of 0.5. Eliminating DLAs with CNR > 3 diminishes the dependence of S on S only slightly, as shown by the red points in the top panel. This is consistent with the small decrease of S for CNR < 3 we saw in Figure 8, which cannot explain the much larger effect of S on the mean S. In other words, the degree of catalogue impurity caused by spectral noise is relatively small, and the large variation of S with S must be caused by a different systematic variation of the selected absorption systems with S . The green points in the top panel show the visually identified DLAs. This visual identification Pâris et al. (2012) has a larger impact reducing the dependence of S on S . The reason is probably that the visual identification was influenced by the presence of metal lines: if metal lines are visually identified, an absorber is more likely to be flagged as DLA.

The middle panel of Figure 10 shows the change in the dependence of the mean S with S when we restrict the sample to DLAs with N HI > 10 20.3 cm -2 (green), and to z < 3.2 (red). As expected, the mean S is increased for higher column density and lower redshift DLAs. There is practically no change in the fractional amount by which the mean S drops with S when eliminating the z < 3.2 DLAs, and a more appreciable change when the low column density DLAs are eliminated. However, the strong decline of the mean S with S is similarly present in all these subsets. The bottom panel compares the corrected S HI dependence on SHI (green), with the previous standard uncorrected case. The correction tends to slightly decrease this dependence.

As we now show, the principal cause of the decline of the mean S with S is that the DLA redshifts are subject to large errors when no metal lines are identified. The Lyα profile of the DLA is often contaminated by neighboring Lyα forest absorption, which may shift the best fit of the DLA redshift. Metal lines, which are narrow, usually pin down the redshift much more accurately. When there is a large redshift error, the metal lines may be partly shifted outside our measurement window, reducing the equivalent width that we compute. This is shown in Figure 11, where the mean S dependence on S is plotted for four widths of the measurement window: 10, 15, 20 and 25 pixels. The drop in the mean S with S is gradually reduced with the measurement window width. For a 25 pixel width, this drop is reduced to only ∼ 10 % for S < 0.5. The redshift error required to shift a metal line by half of the 25 pixel width is ∼ 850 km s -1 , so from the difference in the results for the 20 and 25 pixels, we conclude that redshift errors can reach this value for ∼ 10% of the DLAs. The remaining dependence can be explained by a small fraction of even larger redshift errors and the increase of the mean S with CNR seen in Figure 8.

Catalogue of metal lines and DLA metal strength

We have produced a new catalogue of the DLAs in DR12-DLA, including the metal line equivalent widths computed with the method described in §3. An example of the entries in this catalogue is presented for three DLAs in Table 4: we list the quasar spectrum identification and quasar redshift, the CNR, the DLA redshift and column density, the derived value of S and its error S , the same values corrected for N HI as explained in §4.1, and the equivalent width and errors of each of the 17 metal lines. When lines are not measured for the reasons described in §3.3, the entries for W and its error are set to zero. We make the catalogue publicly available (with metal equivalent widths or not) can be found on the GitHub repository https://github.com/andreuandreu/ Catalogue_Metal_Strength_DLAs_SDSS-BOSS_DR12.

CONCLUSIONS

We have presented a way to classify DLAs according to the strength of their metal lines. A total of 17 metal lines are used to obtain a metal strength parameter that is optimized to be measured with the best possible accuracy. This allows a DLA classification from spectra in surveys with large numbers of objects, where the spectral signal-to-noise is often poor. A catalogue of this parameter and several metal equivalent widths in the DR12-DLA catalogue previously constructed from the BOSS survey using the method of Noterdaeme et al. (2009) is made publicly available.

Our main motivation to present this catalogue is to allow for studies of the dependence of any property of DLAs as a function of the metal strength S, which after correcting for its dependence on N HI to obtain the new parameter S HI , should be a function of the metal abundance and velocity dispersion of the absorbing gas. We plan to do two studies along these lines: first, measuring the bias factor with the technique of Font-Ribera et al. (2012);Pérez-Ràfols et al. (2018) applied to subsamples of the DLAs with different values of S or S HI , and measuring the mean stacked absorption spectrum of the DLA metal lines with the technique of Mas-Ribas et al. (2017) for different values of S or S HI .

Our catalogue presents individual equivalent width measurements, allowing the user to build other combinations of metal line strengths. One of the applications should be to investigate correlations among the strength of equivalent widths of different metal lines. Although the noise is generally very large for individual DLAs, the large size of the sample can allow for detailed studies.

An important limitation in using this catalogue is its imperfect purity. Some of the DLAs may arise due to spectral noise, and some others may be absorbers of low column density with an absorption profile that arises from a velocity distribution of absorbing gas in the Lyα forest but is consistent with a DLA. This can introduce spurious variations of any property we wish to measure as a function of S, if the catalogue purity varies with S. The fact that the mean S as a function of CNR is constant within 5% at CN R > 3 suggests that spectral noise does not introduce substantial purity variations at CN R > 3. We measure a dependence of S on S which we find is mostly due to large redshift errors that are not corrected when no metal lines are significantly detected, which cause the metal lines to be partially shifted out of our measurement window. The impact this

  illustrates the procedure in an example of an Al II 1670 absorption line. The red vertical bars show the limits of the measurement window, and the blue bars are the outer limits of the continuum windows. Green points show the values of the flux in the left panel. The unweighted linear regression is the magenta line. The right panel shows the result of dividing the flux values by this linear regression, giving the transmitted fraction F i . The error on this transmitted fraction, σ Fi , is computed by dividing σ i by the same linear regression.

Figure 1 .

 1 Figure 1. The measurement window (shown by red vertical bars) and two continuum windows (outer limits shown by blue vertical bars) are indicated in this illustrating example of the Al II 1670 line of the DLA in the quasar spectrum MJD-plate-fiber 55182-3587-0100 in BOSS. The flux values in the continuum windows are fitted by a linear regression shown in magenta in the left panel. The ratio of the flux to this linear regression gives the transmitted fraction F i in the right panel.

Figure 2 .Figure 3 .

 23 Figure 2. Histogram comparing the stacked equivalent width from Mas-Ribas et al. (2017) with our mean equivalent width for the 17 lines used in this work, for three different sizes of the measurement window: 10, 15 and 20 pixels. Error bars show the standard deviation of the mean computed from the measured scatter among DLAs.

Figure 5 .

 5 Figure 5. Top panel: Metal strength distribution for three different widths of the measurement window (10, 15 and 20 pixels), compared also with the distribution of S HI for our standard 15 pixel window. Only DLAs with error S < 0.5 are included. Bottom panel: Distribution of the error S for the same three measurement window widths, and the N HI -corrected version for the 15 pixel window.

Figure 6 .Figure 7 .

 67 Figure 6. Green points show the mean metal strength at different column density intervals indicated by the error bars. Blue points show the mean of the corrected metal strength S HI , as defined in 4.1.

Figure 8 .

 8 Figure 8. Top panel: Mean metal strength S as a function of the continuumto-noise ratio in the Lyα forest region. Bottom panel: Mean error S as a function of continuum-to-noise ratio.

Figure 9 .

 9 Figure 9. Upper panel: Distribution of the DLA metal strength for 3 intervals of the error S . Lower panel: Same for the corrected metal strength S HI .

Figure 11 .

 11 Figure 11. Mean metal strength as a function of its error, for four different values of the width of the measurement window: 10, 15, 20 and 25 pixels.

Table 1 .

 1 Metal lines and their principal characteristics used for this work. First column: Name of the metal line. Second column: Wavelength in Angstroms (from Table2 of

Table 1 .

 1 In most cases, there are only 2 or 3 of these highly contributing lines.

		1.0
			All lines
			8 Strongest
		0.8
	Frequency	0.4 0.6
		0.2
		0 0.0	2 Number of W used to compute S 4 6 8 10 12 14	16
	Figure 4. Distribution of the number of lines measured for a DLA among
	all 17 lines in Table 1 (gray solid line), and among the 8 lines with a mean
	contribution C k > 0.15 (defined below in eq. 5).

Table 3 .

 3 Values for the linear regression fitted to the dependence of W on log(N HI ) in equation 6. The first two columns use values of W from the stacking technique shown inMas-Ribas et al. (2017), and the latter 2 columns use W measured with our method, which are used for the analysis in this paper.

	Name	a k stacking	b k stacking	a k	b k
	SiII-FeII 1260	0.316	0.439	0.267	0.568
	OI-SiII 1303	0.591	0.482	0.651	0.489
	CII 1334	0.311	0.443	0.391	0.403
	SiII 1526	0.327	0.270	0.324	0.231
	FeII 1608	0.263	0.089	0.287	0.0744
	AlII 1670	0.338	0.273	0.368	0.233
	SiII 1808	0.101	0.0041	0.110	0.00105
	AlIII 1854	0.091	0.068	0.147	0.0478
	AlIII 1862	0.052	0.038	0.0399	0.0759
	FeII 2344	0.471	0.263	0.568	0.123
	FeII 2374	0.374	0.085	0.337	0.104
	FeII 2383	0.559	0.370	0.513	0.360
	FeII 2587	0.531	0.179	0.574	0.122
	FeII 2600	0.491	0.455	0.528	0.304
	MgII 2796	0.595	0.841	1.008	0.583
	MgII 2803	0.650	0.731	0.923	0.576
	MgI 2853	0.217	0.092	0.208	0.135

dispersion (which affects the degree of saturation of the line), but not on N HI .

  The evolution of the mean metal strength with redshift is shown

		1.20						S HI 15p width
								S 15p width
		1.15
	mean S	1.05 1.10
		1.00
		0.95	2	4	6	8	10 CNR	12	14	16
								ǫ S 15p width
		0.30						ǫ SHI 15p width
		0.25
	ǫ S
		0.20
		0.15
				5			10 CNR		15

  Mean metal strength as a function of the error S for several sub-samples. Upper panel: standard S < 0.5 sample (blue), compared to the sub-sample of DLAs with CNR > 3 (red) and to the visually identified DLAs (green). Middle panel: standard sample compared to DLAs with N HI > 10 20.3 cm -2 (green) and to DLAs with z < 3.2 (red).

			1.4				S 15p width
			1.3				S visual DLA 15p width
			1.2				S 15p width CNR > 3.0
	mean S	0.9 1.0 1.1
			0.8
			0.7
			0.0	0.2	0.4	ǫ S	0.6	0.8	1.0
			1.2 1.3 1.4				S 15p width S 15p width Hd > 20.3 S 15p width z < 3.2
	mean S	0.9 1.0 1.1
			0.8
			0.7
			0.0	0.2	0.4	ǫ S	0.6	0.8	1.0
		1.4				S 15p width
		1.3				S HI 15p width
		1.2
	mean S, S HI	0.9 1.0 1.1
		0.8
		0.7
			0.0	0.2	0.4 ǫ S , ǫ SHI 0.6	0.8	1.0
	Figure 10. Lower panel:
	standard sample compared both for the uncorrected (blue) and corrected
	(green) metal strength.

MNRAS 000, 000-000 (0000)

MNRAS 000, 000-000 (0000)

ACKNOWLEDGEMENTS

We would like to thank George Beker for invaluable guidance and help and Andreu Font-Ribera and Mat Pieri for useful discussions . This work has been supported in part by Spanish grants AYA-2012-33938 and AYA 2015-71091-P. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/.

SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, University of Cambridge, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.

Table 4. Sample of the catalogue for 3 DLAs, with some the information presented: first column redshift of the Quasar, second column Julian(day)-Fiver(number)-Plate(number) from the SDSSIII-BOSS numbering. Third column the DLA redshift. Forth column density as measured by Noterdaeme et al. (2009). Fifth column is the visual flag of the DLA . In columns seven to ten are our corrected and non corrected metal strength (S HI , S) and their errors ( SHI and S as measured in our work. Columns eleven to thirty-three of this table present the measured W and its standard deviation for the sixteen of the seventeen lines (Mg I 2853 is eluded here for presentation reasons). When 0 is expressed for the W and the line has not been selected due to one of the reasons described in §3.3. Other spectrum characteristics are included in the public catalogue but are not shown here for brevity, these are: Balnicity index, Right Ascension, and Declination, following the nomenclature from Noterdaeme et al. (2009)