Detailed 3D structure of OrionA in dust with Gaia DR2

The unprecedented astrometry from Gaia DR2 provides us with an opportunity to study in detail molecular clouds in the solar neighbourhood. Extracting the wealth of information in these data remains a challenge, however. We have further improved our Gaussian Processes-based, three-dimensional dust mapping technique to allow us to study molecular clouds in more detail. These improvements include a significantly better scaling of the computational cost with the number of stars, and taking into account distance uncertainties to individual stars. Using Gaia DR2 astrometry together with 2MASS and WISE photometry for 30 000 stars, we infer the distribution of dust out to 600 pc in the direction of the Orion A molecular cloud. We identify a bubble-like structure in front of Orion A, centred at a distance of about 350 pc from the Sun. The main Orion A structure is visible at slightly larger distances, and we clearly see a tail extending over 100 pc that is curved and slightly inclined to the line-of-sight. The location of our foreground structure coincides with 5-10 Myr old stellar populations, suggesting a star formation episode that predates that of the Orion Nebula Cluster itself. We identify also the main structure of the Orion B molecular cloud, and in addition discover a background component to this at a distance of about 460 pc from the Sun. Finally, we associate our dust components at different distances with the plane-of-the-sky magnetic field orientation as mapped by Planck. This provides valuable information for modelling the magnetic field in 3D around star forming regions.


Introduction
The Orion molecular complex is the nearest site actively forming massive stars in the Galaxy Menten et al. (2007); Bally (2008). As a nearby laboratory, Orion has been studied in various aspects, ranging from distance estimates to different parts of the cloud to the star formation processes and individual stellar populations (e.g. Brown et al. 1994;Menten et al. 2007;Jeffries 2007;Bally 2008;Alves & Bouy 2012;Bouy et al. 2014;Schlafly et al. 2015;Zari et al. 2017;Kounkel et al. 2018;Zari et al. 2019). Stars in the Orion region are also known to be responsible for the creation of the Orion-Eridanus superbubble, a large cavity in the vicinity of the Orion that extends to the constellation of Eridanus in the sky (e.g. Bally 2008;Pon 2015). Bubble structures are very common in the interstellar medium (ISM). They are results of the presence and evolution of young, massive stars that influence their surrounding ISM through radiations, stellar winds and supernovae explosions (e.g. Heiles 1979;Mac Low et al. 1989).
Although the projected picture of the Orion region in the plane of the sky has been determined from various observations of the gas/dust emission (e.g. Ochsendorf et al. 2015;Soler et al. 2018), the distance to different parts of the cloud is still debated. Orion A, the giant molecular filament situated in projection inside Barnard's loop, is possibly the most studied target in this vicinity. Home to the Orion Nebulae Cluster (ONC) at a distance of ∼ 400 pc, Orion A consists of rich clusters of young stars and active sites of massive star formation. Schlafly et al. (2015) demonstrated that Orion A is not a flat filament in the plane of sky. By estimating the distance gradient along the filament, they suggested that the southern part of the filament (hereafter tail) is further away from the Sun than the northern part (hereafter head) hosting the ONC. Recent work by Großschedl et al. (2018) revealed an extended tail of Orion A to larger distances than previously estimated based on the distribution of the young stellar objects (YSOs) in the Orion A vicinity. Zucker et al. (2020) confirmed the distance gradient for Orion A by estimating the distance for individual sight-lines along the filament.
One way of probing the structure of the cloud is through mapping its full three-dimensional (3D) dust distribution. Dust, only a tiny fraction of matter in the Galaxy, scatters, absorbs and re-emits light making distant objects look fainter than they are. In addition to having negative effects on observations of more distant objects, dust plays important roles in creating and shaping the ISM. It protects molecules from the high energy UV radiation, which would otherwise impede star formation. The ISM and the life cycle of stars are tightly related. Therefore studying different properties of the ISM, including the 3D distribution of dust in the Galaxy, can provide valuable information about the structure of the ISM and potential sites of star formation in the Milky Way. Furthermore, the 3D distribution of dust provides crucial information to model the distribution and dynamics of star-forming clouds. For example, knowledge of distances to various dust overdensities along a line of sight (l.o.s) is helpful to disentangle the components responsible for the plane-of-thesky magnetic field orientations derived from submillimeter polarization observations, such as those by ESA's Planck satellite (Planck Collaboration Int. XXI 2015).
Various dust extinction mapping techniques have been developed over the past few years. While most of these methods infer the individual line-of-sight extinction towards stars, then take the derivative of these individual sight lines to get the extinction per unit distance (e.g. Marshall et al. 2006 2018), dramatic artefacts produced in these approaches make it hard to explore the physical properties of the ISM, like the structure of the molecular clouds. Green et al. (2019) tried to overcome this drawback by applying a smoothing function on the extinction derivatives using an iterative approach. Lallement et al. (2019) focused on mapping the differential extinction in 3D by taking into account the neighbouring correlations. Although this approach has the advantage of producing smooth maps, it does not consider distance and extinction uncertainties when inferring the differential extinction, therefore could be biased due to the data quality cuts.
Probing the 3D distribution of the dust towards the Orion region can provide valuable information about the distance to and structure of the cloud. The 3D dust mapping towards the Orion complex by Schlafly et al. (2015) revealed the Orion dust ring. To further study the dust distribution towards Orion, as shown in (Rezaei Kh. et al. 2017, 2018b, we have developed a nonparametric 3D dust mapping technique that takes into account the 3D correlation between points in space, allowing the model to trace arbitrary dust variation. In (Rezaei Kh. et al. 2018a) we mapped the dust distribution towards the Orion region and demonstrated the capability of our method to capture different dust clouds in the 3D space without having the discontinuities and artefacts presented in other works. The latest data release from the second release of the Gaia satellite (Gaia DR2, Gaia Collaboration et al. 2018) with unprecedented astrometry enables us to further study this area.
In the present work we improve our mapping technique by including both distance and extinction uncertainties, together with resolving the computational limitations, thereby allowing us to exploit a large dataset like Gaia DR2 as the input, and to produce a detailed 3D dust map of the Orion A. We also compare the location of our dust clouds with that of young stellar populations in the region, then present an analysis of the magnetic field orientation using our density structures.

3D mapping technique
Here we briefly summarise our modelling approach (refer to Rezaei Kh. et al. (2017) and Rezaei Kh. et al. (2018b)) and also explain the improvements made.
We use a non-parametric model to infer the local dust density using the position and l.o.s attenuation to a number of stars in 3D space. The attenuation (a) is related to the extinction (A) as A 1.0857a. We divide each l.o.s into small 1D cells to approximate the attenuations as the sum of the dust densities in the cells along the l.o.s to the stars. Afterwards, we connect all these cells in the 3D space using a Gaussian process that takes into consideration the neighbouring correlations between points. The closer two points are in the physical space, the more correlated their local dust densities. We use a truncated covariance function in order to account for the correlations, i.e. the points will be correlated only if they are closer than a correlation length. This way we set up our model to then infer the dust density for any arbitrary point in this space, even along a l.o.s that was not initially observed. In addition to the specification of the cell size (with fixed length), our model has three hyper-parameters: λ which is the correlation length, θ which sets the amplitude of the density variance, and the mean of the Gaussian process. We set these hyper-parameters based on the input data (using the approach described in Rezaei Kh. et al. 2017, 2018b.
As we explained in the aforementioned papers, there were a couple of limitations in our model that we overcome in the present work. One of the main drawbacks of the model was that we did not take into account distance uncertainty. This limited us to only using data with very precise distance measurements. Since the position of input stars and predicting points are "given" in our model (Rezaei Kh. et al. 2017, 2018b, any direct application of the distance uncertainty in the analytic solution is not possible. We have overcome this issue by propagating the distance uncertainty of a star into its attenuation uncertainty. This way we have where a is the attenuation, ρ is the mean dust density along the l.o.s, σ r is the distance uncertainty, and σ a d is the uncertainty in attenuation propagated from the distance uncertainty. The total input attenuation uncertainty in the model (σ a tot ) is then where σ a is the measured attenuation uncertainty. This provides us with an upper limit on the input attenuation uncertainty as a result of the distance uncertainty: if a star is in or close to a high density environment, then varying its distance would cause a significant difference in the attenuation measurements; otherwise, if changing the star's location would not impact its attenuation measurement, this would overestimate its attenuation uncertainty. We use this maximum possible uncertainty as the input to our model to account for both the extinction and distance uncertainties. This provides, in the higher-density regimes, more flexibility for the model to capture the underlying dust density variations when considering the neighbouring correlations.
Another major improvement we report in this paper is on the computational limitation. This was an issue we discussed in Rezaei Kh. et al. (2017) and Rezaei Kh. et al. (2018b): since we consider the correlation between all points in space, the further the stars, the more cells we have, resulting in a more computationally expensive calculation. We have overcome this problem by partitioning our dataset into different slices as we now explain: First, we infer the dust density for an inner sphere with radius of one correlation length, λ, using stars within two correlation lengths. We then treat the inferred inner densities as known values, taking into account their correlated inferred density uncertainties, to derive the densities for the next layer and continue until we reach our desired distance. From the second layer on, points in different directions could be more distant than λ from the l.o.s. So to predict dust densities along a given l.o.s, we only consider stars that lie within the correlation length of that l.o.s. This corresponds to a cylinder with radius equal to λ (extending from the sphere of radius λ that was the first layer). This approach ensures that we include all relevant neighbouring stars for the correlation computations, yet speeds these by ignoring stars beyond a distance λ.
Using this setup we can predict the dust densities for any point in space. The resolution of the final map is set by the density of the input data, which is on average equal to the typical separation between the input stars.

Input Data
We use data from the second Gaia data release (Gaia DR2, Gaia Collaboration et al. 2018) to get the 3D positions of stars. It is important to note that Gaia provides parallaxes for stars and not The predictions are made on regular grids for every 0.5 degrees in the Galactic l and b, and every 10 pc in distance. The 2D image is then produced by applying a smoothing kernel (with 4-pc scale length) to handle the missing pixels. In order to not produce extra smoothing than that of the method, the length scale of the kernel is chosen to be relatively small; hence, the distance gridding is still apparent in the left panel.
distances and since we do not cut on the precise parallax measurements, we need to account for the noisy parallaxes to get the distance information. We therefore use the catalogue of Bailer-Jones et al. (2018) who infer distances to Gaia sources from noisy parallaxes. For the uncertainty in the estimated distance, we take the average of the estimated lower and upper confidence interval.
In addition to 3D positions, the method of course needs a measure of extinction towards each star. Gaia DR2 provides extinction measurements for around 88 million sources (Andrae et al. 2018). However, due to the strong degeneracy between the extinctions and temperatures of stars, the individual extinctions from Gaia DR2 are less reliable for star-by-star analyses . We, therefore estimate extinctions using the Rayleigh-Jeans Colour Excess (RJCE, Majewski et al. 2011) method, which uses near-infrared minus midinfrared (NIR-MIR) colours of stars in order to get their K s -band extinctions. The RJCE method works based on the fact that the distribution of the intrinsic colours of stars in NIR-MIR is so narrow that it can be treated as a known value with some uncertainty. The difference between the observed colours and this intrinsic value is then due to the column of the dust between the star and the observer (Majewski et al. 2011). We use H band photometry from the Two Micron All-Sky Survey (2MASS; Skrutskie et al. 2006) as the near-infrared data and the WISE W2 band photometry (Wright et al. 2010) as the mid-infrared one. Both catalogues are cross-matched with the Gaia DR2 sources in the Gaia archive 1 .
In this work, we focus on the Orion A region and we select stars within 204 • < l < 218 • and −22 • < b < −15 • . As explained in detail in Rezaei Kh. et al. (2018a), after we derive the extinction values using the RJCE method, we select our final sample based on the positions of the stars on the de-reddened colour-magnitude diagram in order to remove the outliers that happen to get unrealistic high extinction values due to the star being photometrically variable, or so young that it still harbours a dusty disk (see Rezaei Kh. et al. 2018a). This leaves us with Fig. 1 shows 2D Cartesian projections of the inferred 3D dust densities in the Orion A region. As can be seen from the figure, there is a bubble-like structure in the foreground of the main Orion A filament, at about 350 pc, that expands to larger longitude and latitude compared to Orion A. In addition, the extent of the Orion A to further distances is clear from the left panel. The head of the Orion A, which includes the Orion Nebula Cluster (ONC), appears to be at around 400 pc, while as seen from both projections, the tail of the cloud is extending to distances of about 490 pc, making the total length of the Orion A cloud to be over 100 pc.

Orion dust map
The distance to and the location of individual parts are better seen in Fig. 2 that shows the observer's view of the same region. Each panel indicates a slice through the region at a fixed distance. The foreground structure is clearly seen in the 345-pc panel and seems to be extended upwards to around l = −16 • . The extent of the tail of the Orion A is also visible in multiple panels up to further distances.
Our map covers a slightly larger area than that of the Orion A; consequently, the lower part of the Orion B (higher densities above latitude of −17 • ) appears in our results. As shown in figures 1 and 2, apart from the main Orion B over-density at slightly above 400 pc, a background component is revealed at the distance of around 460 pc. A dedicated study on the Orion B cloud will be carried out in a future work.
In the previous figures we plotted only the mean of our density predictions, while our method provides the probability density distribution for each point. To investigate the uncertainties of the predictions and the significance of the predicted dust clouds we look at two lines-of-sight towards the foreground cloud and plot the predicted dust densities and their uncertainties as a function of distance, as illustrated in Fig. 3. The two l.o.s are also shown on the 3D projections in Fig. 1 Fig. 3. Dust density vs. distance for two different l.o.s. towards upper (blue shades) and lower (grey shades) parts of the foreground structure (overplotted on Fig. 1). The black line shows the mean and the shades represent one standard deviation (also computed by the Gaussian Process model). The grey-shaded first peak represents the foreground cloud while it seems to be connected to the main cloud in the background (the second peak). The two blue peaks illustrate the front and back edges of the foreground cloud.
foreground cloud is evident in both lines-of-sight as significant overdensities. As demonstrated in the plot, the lower part of the foreground structure (grey) is connected to the main Orion A filament at larger distances, while the density of the upper part (blue) drops significantly after the foreground cloud, indicating the separation between the two at the upper part.
It is important to note that our model infers densities for "given" points in space, i.e. the location of the predicting points are fixed. As a result, the model by design does not provide distance uncertainty for the predictions. The model does, however, provide uncertainties for the inferred densities. Similar to the way we accounted for the distance uncertainty in the input data (see eq. 1), inferred density uncertainties can be translated into distance uncertainties. This way the typical distance uncertainty is estimated to be 10 pc in our map.

The Orion A structure and young stellar populations
The foreground bubble in our map has not been reported prior to this work. This could be mainly due to the projection effects, lack of distance resolution, or the limitations of the underlying techniques. From about 30,000 stars in our sample, nearly 5,000 of them are within 300 -400 pc distance, of which ∼ 2200 are located around the coordinates where we discovered the foreground cloud. This indicates the high number statistics involved in our analysis.
To elaborate this more, we compare the recent results of Großschedl et al. (2018) Fig. 4. Same as Fig. 3 but for three other l.o.s; one towards the Orion A head (in grey), another one towards the middle of the filament, at the location of the L1641-S cluster (in blue), and the third one towards the tail of the Orion A and at the location of the L1647-S cluster (in red). The projected sight lines are overplotted on Fig. 2 as a reference. the average distances of stars in 1-degree longitude bins. They have reported that the Orion A cloud is not a straight filament as is seen in 2D projection, but has an elongated tail extending to larger distances (Großschedl et al. 2018). The tail of Orion A is clearly seen in the left panel of Fig. 1 in our work, here we focus on three specific l.o.s as probed by Großschedl et al. (2018) which we show in Fig. 4: one towards the head of the Orion A, and the other two towards the two star clusters L1641-S and L1647-S. Going from the head to the L1641-S and L1647-S, the peaks of the density distribution moves from smaller to larger distances, confirming the results of Großschedl et al. (2018). We note that the l.o.s towards the head of the Orion A shows two peaks, the closer one being the foreground cloud (as was discussed in the previous figures/sections). We also see that the dust density towards the tail of the cloud shows a wide distribution meaning that the tail starts as close as about 400 pc, then extends to further distances. These points are indeed present in Fig. 3 of Großschedl et al. (2018): at longitudes of around 208 • to 210 • , there are quite a few YSO distances around 350 pc, where we discovered the foreground cloud. In addition, the large scatter in their distance estimates around the middle and tail of the filament is indeed real representing the extent of the cloud. Fig. 5 shows the YSOs in Großschedl et al. (2018), over-plotted on our distance slice plot. The YSO locations nicely match our dust distribution in the regions covered by Großschedl et al. (2018); both at the location of the foreground cloud, and along the main Orion A filament.
The location of our foreground structure seems to be consistent with the findings of recent studies focused on the stellar populations towards Orion A. Bouy et al. (2014) reported a large foreground population towards Orion A, loosely clustered around NGC1980 and NGC1981, at a distance of 380 pc. They estimated the age of this population to be between 5 and 10 Myr. Kounkel et al. (2018) reported similar results for their Orion D structure where, especially around NGC1981, stars are found at a distance of 375 ±3 pc. Jerabkova et al. (2019) also confirmed the presence of multiple stellar populations towards the ONC and indicated that the older population is closer to us than the ONC. Zari et al. (2019) confirmed that the Orion OB association consists of numerous groups with different ages and kinematics and suggested a complex star formation history. In particular, their B1 group (distance of 365 ±10 pc, age of 10±0.5 Myr) appears to be located around our reported foreground dust structure. The fact that the location of our foreground dust cloud and that of older stellar populations roughly coincide, suggests that star formation in the Orion region might indeed have started at a closer distance to the Sun than that of the ONC. Supernova explosions and stellar winds from massive stars could have triggered a new episode of star formation in the ONC and other younger clusters in the region. Our foreground dust cloud could then be a remnant of a previous star formation episode.

Reconstructed extinction map
We can test our inferred 3D distribution of the dust in this region by recreating a 2D extinction map. We do so by summing up the predicted densities along each l.o.s. from the Sun to 600 pc, as presented in Fig. 6. We superimpose the borders of Orion A and Orion B from the Lombardi et al. (2011) extinction map. The main Orion A structure can be seen in the central part, together with the lower part of the Orion B that appears at the top of the plot. Around the longitude of 214 • , a ring-type structure is visible. This is also seen in the dust column density map from Planck (see Fig. 7). The ring was thought to be associated with the κ-Ori star located at the center of the ring in the 2D projection. The distance to the κ-Ori star is estimated to be 200 ±10 pc from Hipparcos (van Leeuwen 2007); however, there does not seem to be any associated overdensity in our map at those distances. Since the star is fairly bright (V ∼ 2 mag), Gaia has not yet published its parallax. We do not find a ring structure in our 3D map that would correspond to what is seen in the 2D projection. More focused study on the κ-Ori region will be carried out in future work.
It is important to note that we use optical data from Gaia, which omits stars in very dense parts of the cloud due to the high obscuration by dust. When computing the correlations, all stars within a correlation length contribute to the inference of the density for a point in that volume. Since lower extinction stars are represented more in the input data, the average extinction in dense regions is less than it would be in the presence of deeper dataset. Therefore, the values we report here for extinctions in the higher density regime (like around the ONC) could be an underestimate of the total extinction.

Magnetic field structure
One application of our study is the analysis of the plane-ofthe-sky magnetic field orientation (ψ) revealed by the polarized emission from interstellar dust grains at submillimeter wavelengths (see, Andersson et al. 2015, for a review). The magnetic field orientation inferred from the observations by Planck at 353 GHz toward Orion (Planck Collaboration Int. XXI 2015), shown in Fig. 7, reveals a variety of bends and kinks that can be attributed to the combined effect of magnetic tension, gravity, and turbulence (see for example, Planck Collaboration Int. XX 2015; Planck Collaboration Int. XXXV 2016). However, it is also possible that some of the features in the map are the sum of the contributions from the magnetic fields in dust overdensities located at different distances, which are overlapped in projection.
The observed polarization can be modelled as a series of dust slabs located at different distances, each with its own mean magnetic field orientation, like layers in a diorama (Martin 1974). 365 pc q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q qq q q qq q q qq q q q qq q q q q q q q q q q q q q q q q q q q qq q q q q q q q q qq q q q q q q q q qq q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q Given that the magnetic field is a vectorial quantity, the contribution from each slab to the total polarization is not proportional to the amount of dust in each slab but depends on the orientation of the magnetic field. We evaluate here whether the slab model can provide further information on the structure of magnetic field toward Orion A by considering the field orientations within dust overdensities at different distances. To do this, we construct 16 slabs of extinction between 180 and 600 pc by integrating the Fig. 7. Dust column density and magnetic fields toward the Orion clouds. The colors represent the logarithm of the dust column densities inferred from the Planck observations (Planck Collaboration XI 2014). The overlayed drapery pattern correspond to the representation of the plane-of-the-sky magnetic field orientation derived from the Planck polarization observations at 353 GHz obtained using the line integral convolution method (LIC, Cabral & Leedom 1993). The green contours correspond to the areas defined by the extinction threshold A Ks ≥ 0.06 mag in the extinction slabs at the indicated distances. dust densities with ±15 pc of the distances shown in Fig. 5. For each distance slab, we identify the region in which the extinction is above a specified threshold and calculate the mean orientation of the magnetic field ψ within that area. The changes in ψ in different slabs indicates that there are particular dust overdensities that might be responsible for the observed field orientation.
If the values of ψ are roughly constant across distances, this method is insufficient to relate the field orientations to the dust structure in 3D. Figure 8 shows the values of ψ in the regions defined by five different extinction thresholds between A K s = 0.04 and 0.08 mag in the slabs at six different distances. We see that for the extinction thresholds A K s < 0.06 there is no significant variation in ψ , most likely because the magnetic field orientation is being averaged over similarly large areas on the plane of sky for all distance slabs resulting in similar values for ψ . However, for the extinction thresholds A K s ≥ 0.06 we see changes of up to 50 • within the regions defined by the extinction thresholds at different distances. The contours in Fig. 7 illustrate an example of regions where the most significant changes in ψ are found for the particular case A K s ≥ 0.06. The figure clearly shows that the mostly vertical magnetic field orientation at l = 209 • , b = −20 • is associated with the dust overdensities in the slabs centred around 385 pc. In contrast, the mostly horizontal magnetic field at l = 214 • , b = −19 • corresponds to the overdensity in the slab corresponding to 405 pc, as illustrated in Fig. 9. The region of overlap between the 385-pc and the 405-pc thresholds shows a field orientation at around 45 • . This could be the result of the superpositions of the orientations in the two slabs, although the current analysis cannot provide further evidence of that scenario.
The fact that different dust components along a l.o.s are causing changes to the magnetic field orientation is important for three reasons. First, it indicates that the models of the magnetic field orientation, such as those presented in Heiles (1997) and Tahani et al. (2019), are limited when considering Orion A as just one object. They would instead need to account for the fact that the field is changing its orientation in a dust structure that is just 30 pc behind the main body of the cloud. Second, it indicates that the field is bending around and within the cloud. This suggests that this structure is not fully dominated by one single 10-pc-scale magnetic field but rather is showing the combined effects of trans-Alfvénic turbulence (i.e. kinetic and magnetic energy densities are comparable) and gravitational collapse (Hennebelle & Inutsuka 2019;Pattle & Fissel 2019). Third, it provides an example of a technique that can complement and Fig. 9. Extinction in the 30-pc slabs around the distances 385 and 405 pc, shown in red and cyan as a two-colour image. Overlap in the extinction in the two slabs is show by the brighter colour between cyan and red (e.g. around l = 213 • , b = −19 • ). The contours correspond to the logarithm of the dust column densities inferred from the Planck observations, log 10 (N H /cm) −2 . The drapery pattern corresponds to the orientation of the plane-of-the-sky magnetic field orientation inferred from the Planck 353-GHz polarization observations. enhance the studies of the magnetic field in 3D using polarized starlight observations (see for example Panopoulou et al. 2019).

Concluding remarks
We have developed a sophisticated mapping technique that takes into account the neighbouring correlation in the space using a Gaussian process prior with a truncated covariance function. The method also considers the distance and extinction uncertainties to individual stars, enabling us to provide a detailed 3D dust map of local molecular clouds. We have presented the results towards the Orion A region where, for the first time, we reported a bubble-like dust overdensity at 350 pc in front of the previously known Orion A cloud. We also illustrated the whole shape of the cloud with its extended tail to distances of about 490 pc. This indicates a length of more than 100 pc for the filament. We also reported a background component to Orion B at a distance of 460 pc.
The presence of stellar associations older than that of the ONC around the same location as our foreground overdensity suggests an early star formation episode in front of the ONC. This could have triggered subsequent episodes of star formation in the region.
We also analysed the magnetic field orientation in the plane of the sky. We connected variations in the magnetic field orientation angles with variations in dust density along the l.o.s. This can provide valuable information for understanding the magnetic field distribution in 3D around a star-forming region, which is a crucial step forward in our understanding of the role of the magnetic field in the process of star formation.