Observationally derived magnetic field strength and 3D components in the HD 142527 disk

HD 142527 is a Herbig Fe star¹ with a mass and age of ~2.4 M_⊙ and ~3 Myr, respectively^2,3,4, and it is 157 pc away from our Solar System^5,6. The accretion rate to the star has been estimated to be 10⁻⁷ M_⊙ yr⁻¹ (ref. ⁷). The disk surrounding the young star exhibits a wide gap with a radial width of about 100 au and a horseshoe-shaped ring structure^8,9. The high-contrast emission from its north side is brighter than from its south side, which suggests that there are large dust grains in the northern region and small dust grains in the southern region^10,11. Therefore, the protoplanetary disk around HD 142527 is an ideal laboratory for investigating grain growth, which is considered to be the first step in planet formation. We analysed the polarization data of the dust continuum emission at wavelengths of 0.87 mm (band 7), 1.3 mm (band 6), 2.1 mm (band 4) and 2.7 mm (band 3) obtained by the Atacama Large Millimeter/submillimeter Array (ALMA) (Methods).

Figure 1 shows images of Stokes I, polarized intensity (PI) and polarization fraction (PF) at these wavelengths. The polarization vectors (PA) are also plotted on the PI and PF images. The horseshoe structure can be recognized in the Stokes I images at all wavelengths, as revealed in previous observations^8,9,10,11. Although the Stokes I images are similar in these wavelengths, the polarization images show different morphologies. In particular, the polarization vectors and fraction in the northern region change significantly with wavelength. The 0.87 mm polarization vectors show a flip around the peak of the Stokes I emission, as also shown in previous studies^10,12, whereas the 2.1 and 2.7 mm polarization vectors seem to be almost azimuthally aligned. The 1.3 mm polarization vectors also show the flip, although it is less prominent¹³. These different polarization patterns can be explained by the contribution of the self-scattering, which is efficient only when the maximum grain size is close to the observing wavelength. Therefore, dust grains of size ~100 μm are abundant in the northern region^12,14. In contrast to the northern region, the southern region shows almost the same polarization pattern regardless of the observed wavelength. All the polarization vectors are in the radial direction, and the polarization fraction reaches as high as 15%.

Stokes I (top), polarized intensity (PI) (middle) and polarization fraction (PF) (bottom) at 0.87, 1.3, 2.1 and 2.7 mm wavelengths plotted against right ascension (Δα) and declination (Δδ). Polarization vectors are also shown on the PI and PF images. The Stokes I and PI are converted to the brightness temperature (T_B) by assuming the Rayleigh–Jeans approximation. The PF and polarization vectors are plotted where the polarized intensity is detected above the 3σ noise level. Contours on the PI and PF images indicate the brightness temperatures of the Stokes I emission of [0.1, 1, 10, 20] K for 0.87 mm, [0.1, 1, 5, 10, 20] K for 1.3 mm, [0.1, 1, 5, 9] K for 2.1 mm and [0.06, 0.3, 1, 3] K for 2.7 mm. The resolution (beam size) is shown as a small ellipse in the lower left corner of each panel.

Full size image

The consistent polarization vectors and polarization fractions in the southern region (position angle of 100–250°; Methods and Supplementary Figs. 1, 2 and 3) suggest that the polarization is produced by grain alignment rather than self-scattering. The high polarization (~15%) in the radial direction is expected to originate from the magnetically aligned grains^15,16 under the toroidal magnetic field^17,18,19. By contrast, other alignment mechanisms, such as radiative alignment or mechanical alignment, predict azimuthally polarized emission^20,21,22, which may be the case for the 1.3, 2.1 and 2.7 mm polarizations in the northern region. Note that the magnetic alignment also predicts azimuthally polarized emission at millimetre wavelengths when the dust is larger than 1 mm due to the negative polarization fraction in the Mie regime²³. Therefore, we conclude that at least the smaller dust grains in the southern region are magnetically aligned, as occurs in the interstellar medium^15,16. This may imply that the dust grains in the southern region contain large iron inclusions that are aligned with the magnetic field in the protoplanetary disks²⁴.

For the southern region, the polarization vector is expected to be perpendicular to the magnetic field direction projected onto the sky plane. Thus, we investigated the magnetic morphology by rotating the polarization vectors by 90°. Figure 2 indicates that the toroidal field is dominant, as expected from magnetohydrodynamics simulations^25,26,27. To further constrain the magnetic field structure, we measured the angular offset of the observed magnetic field from the purely toroidal case, assuming a disk inclination and position angle of 27° and 161°, respectively^6,28.

The magnetic field is shown by the white segments superimposed on the 0.87 mm Stokes I continuum image, which was derived by rotating the 0.87 mm polarization vectors by 90° because the dust polarization is caused by the thermal emission of magnetically aligned dust grains in the southern region. The black dashed line indicates the azimuthal direction of the disk and has a radius of 210 au.

Full size image

Figure 3 shows the angular difference (Δϕ = B_ϕ − B_obs) between the observed magnetic field (B_obs) and the pure toroidal magnetic field (B_ϕ) as a function of position angle at 200 ± 50 au from the central star in the southern region. From this figure, we found that the observed magnetic field has systemic offsets of about 10° from the toroidal case. We interpret these offsets as being caused by the radial and vertical components of the magnetic fields in the disk.

Angular difference (Δϕ = B_ϕ − B_obs) between the observed magnetic field (B_obs) and the pure toroidal magnetic field (B_ϕ) as a function of position angle at 200 ± 50 au from the central star in the southern region. The colour indicates the radius from the central star. The data were obtained using Nyquist sampling (half of the beam size) with a disk radius of 200 ± 50 au. The data are presented as the median of the distribution with the errors corresponding to the 68th percentiles. There is a systemic offset of about 10°, which was caused by the contributions of the radial (B_r) and vertical (B_z) magnetic fields. The black line shows the best-fitting parameter ∣B_r∣:∣B_ϕ∣:∣B_z∣ = 0.26:1:0.23 that explains this offset. Note that the fitting errors are not shown because they are not meaningful due to the large, reduced chi-square value of approximately 9 for the best-fitting model.

Full size image

Single-layer model

Here we propose a method for estimating the relative strengths of the three-dimensional components of the magnetic fields: radial (B_r), azimuthal (B_ϕ) and vertical (B_z). The B_r and B_z components act to deviate the sky-projected magnetic field vector from the purely toroidal direction. Their effects vary at different locations on the disk (Fig. 4), causing azimuthal variations in the angular offsets, even if the intrinsic magnetic field exerted on the disk is homogeneous. The angular offset was estimated with

where ({({B}_{r})}_{phi perp }) is the component perpendicular to the azimuthal direction in sky-projected B_r, ({({B}_{r})}_{phi }) is the component parallel to the azimuthal direction in sky-projected B_r, ({B}_{z,mathrm{proj}}={B}_{z}sin {rm{(inclination)}}) is the sky-projected strength, ({({B}_{z,mathrm{proj}})}_{phi perp }) is the component perpendicular to the azimuthal direction in B_z,proj and ({({B}_{z,mathrm{proj}})}_{phi }) is the component parallel to the azimuthal direction in B_z,proj. These parameters are functions of the position angle, as can be seen by considering the projected direction of the grain alignment in the inclined disk geometry. Thus, the variation of the offset value Δϕ with position angle allowed us to estimate the relative strengths of the three-dimensional magnetic field. We found that ∣B_r∣:∣B_ϕ∣:∣B_z∣ = 0.26:1:0.23 from the chi-squared test best fitted the observations, as shown in Fig. 3. The reduced chi-square value for the best-fitting model is approximately 9, suggesting that more detailed observations and refined models are necessary for giving a better understanding of the fine magnetic structure, such as the radial distribution. Nevertheless, overall, the absolute values and variation of Δϕ suggest that the azimuthal (toroidal) magnetic field is dominant, whereas the radial and vertical magnetic field strengths are almost comparable. Note that there is uncertainty in the magnetic field direction because the polarization vectors indicate only the alignment direction. Therefore, we used the absolute values for the relative strength.

We propose that the dependence of the angular offsets (Δϕ) on the disk position angle (shown in Fig. 3) is caused by the three-dimensional magnetic field structure of the disk as B_r, B_ϕ and B_z. At position P, the radial and vertical magnetic fields point in the same direction, strengthening the offset from the azimuthal direction, whereas at position Q, the radial and vertical magnetic fields point in opposite directions, weakening the offset from the azimuthal direction. Therefore, Δϕ varies with the disk position angle.

Full size image

We assessed the dependence of the relative strengths ∣B_r∣:∣B_ϕ∣:∣B_z∣ on the angular offsets (Δϕ) by changing these parameters to confirm the robustness of the best-fitting results. Figure 5 shows four different cases for the relative strengths of the magnetic field. The black line indicates the best-fitting ratio, and the dashed coloured lines are for relative strengths ∣B_r∣:∣B_ϕ∣:∣B_z∣ = 0.3:1:0.1 (red), 0.25:1:0.5 (blue) and 0.2:1:0.3 (green). As shown in this figure, the systemic offset of Δϕ is sensitive to B_r, whereas the variation of Δϕ with position angle is sensitive to B_z. Therefore, the absolute values of Δϕ and its dependence on the position angle allowed us to derive the relative strengths of B_r, B_ϕ and B_z. By comparing these different parameters, we suggest that the best-fitting ratio ∣B_r∣:∣B_ϕ∣:∣B_z∣ = 0.26:1:0.23 is a reasonable estimate.

Angular difference (Δϕ = B_ϕ − B_obs) between the observed magnetic field (B_obs) and the pure toroidal magnetic field (B_ϕ) as a function of position angle at 200 ± 50 au from the central star in the southern region. The colour indicates the radius from the central star. The data were obtained using Nyquist sampling (half of the beam size) with a disk radius of 200 ± 50 au. The data are presented as the median of the distribution with the errors corresponding to the 68th percentiles. There is a systemic offset of about 10°, which is caused by the contributions of the radial (B_r) and vertical (B_z) magnetic fields. The black line shows the best-fitting parameter ∣B_r∣:∣B_ϕ∣:∣B_z∣ = 0.26:1:0.23 that explains this offset. The dashed coloured lines indicate the different relative magnetic field strengths, ∣B_r∣:∣B_ϕ∣:∣B_z∣ = 0.3:1:0.1 (red), 0.25:1:0.5 (blue) and 0.2:1:0.3 (green), which shows the dependence of Δϕ on the relative magnetic field strengths.

Full size image

Dual-layer model

It could be expected that B_ϕ and B_r are probably in opposite directions with respect to the disk midplane due to the entrainment of the initial B_z field by the disk rotation (Extended Data Fig. 1). Thus, the observed offset may be due to a combination of an upper disk layer (Δϕ₊) and a lower disk layer (Δϕ₋) (dual-layer model). By considering this effect, we fitted the offset values as a function of the position angle (Methods). The best fit resulted in ∣B_r∣:∣B_ϕ∣:∣B_z∣ = 0.21:1:<0.1, as shown by the black line in Fig. 6. The reduced chi-square value for the best-fitting model is approximately 16, which is larger than that of the single-layer model. This is because the variation of (overline{Delta phi }) (a combination of Δϕ₊ and Δϕ₋; equation (2)) with position angle cannot be reproduced by the dual-layer model. To confirm this, different magnetic field strengths are also plotted (coloured dashed lines in Fig. 6). We found that larger B_z causes a smaller angular offset of (overline{Delta phi }) around the disk major axis (position angle of 150°) and a larger angular offset around the disk minor axis (position angle of 220°). This is because B_z has the same direction as B_ϕ along the disk major axis (which suppresses (overline{Delta phi })) and as B_r along the disk minor axis (which enhances (overline{Delta phi })). We also found that B_z becomes more efficient when it has the same direction as B_ϕ than when it has the opposite direction to B_ϕ in the dual-layer model because the projected length of the elongated dust grains becomes longer. By contrast, the previous single-layer model can enhance Δϕ in the disk major axis and suppress Δϕ in the disk minor axis by changing the direction of B_z with respect to B_r. Therefore, the best-fitting parameter of the single-layer model suggests that the B_r and B_z directions are opposite to each other. That is, if B_r is outward, then B_z is from south to north in the sky-projected direction.

Angular difference (Δϕ = B_ϕ − B_obs) between the observed magnetic field (B_obs) and the pure toroidal magnetic field (B_ϕ) as a function of position angle at 200 ± 50 au from the central star in the southern region. The colour indicates the radius from the central star. The data were obtained using Nyquist sampling (half of the beam size) with a disk radius of 200 ± 50 au. The data are presented as the median of the distribution with the errors corresponding to the 68th percentiles. There is a systemic offset of about 10°, which is caused by the contributions of the radial (B_r) and vertical (B_z) magnetic fields. For the dual-layer model of the magnetic field, the best fit gives ∣B_r∣:∣B_ϕ∣:∣B_z∣ ≈ 0.21:1:<0.1, as shown by the black line. The different relative strengths are plotted as the coloured dashed lines: ∣B_r∣:∣B_ϕ∣:∣B_z∣ = 0.21:1:0.4 (red), 0.21:1:0.6 (blue) and 0.21:1:1 (green).

Full size image

These models imply that the single-layer model can reproduce the observed features of the systemic offset and the variation with position angle better than the dual-layer model. This suggests that, due to turbulence, the magnetic field across the disk midplane is not purely symmetric. Some magnetohydrodynamics simulations have shown such asymmetric magnetic field structures²⁹, although few simulations have been performed for disks with substructures such as rings and crescents³⁰. It is possible that the magnetic field in a crescent disk like this one is not a simple symmetric structure across the midplane. It will be necessary to investigate in further simulations the magnetic field in ring or crescent structures. Note that the emission was estimated to be optically thin (τ ≈ 0.1) in the southern region at a wavelength of 0.87 mm, suggesting that the polarized emission is the average value along the vertical direction of the scale height and mainly traces the midplane region. This also indicates that the single-layer model is better than the dual-layer model, because B_ϕ is assumed to fall at the disk midplane in the dual-layer model so that it can flip direction, whereas the observations indicate that B_ϕ is the dominant component of the disk magnetic field. Although different values of B_z were estimated by the two models, B_r has almost the same strength of B_r ≈ (0.2–0.3)B_ϕ in both. This is because B_r was mainly determined by the systemic offset of Δϕ, whereas B_z was determined by the variation of Δϕ with the position angle. Therefore, B_r could be estimated better than B_z, and B_r ≈ (0.2–0.3)B_ϕ is robust. As the reduced chi-square values for both the single- and dual-layer models are much larger than 1 (9 for the single-layer model and 16 for the dual-layer model), more detailed observations and refined models are necessary to provide a better understanding of the fine magnetic structure, such as the radial distributions of B_r and B_z.

Our results for the relative magnetic field strength allowed us to estimate important physical parameters related to the disk magnetic field based on observations. To realize B_z/B_ϕ ≈ 0.23 in the disk around HD 142527, we estimated the ambipolar resistivity η_A ≈ 3.0 × 10¹⁷ cm² s⁻¹, the magnetic field strength B ≈ 0.3 mG, the ambipolar Elsasser number Am ≈ 0.4 and the plasma beta β ≈ 2.0 × 10² (Methods). The Am parameter is a dimensionless number that characterizes the ambipolar diffusion. From these estimates, we suggest that magnetorotational instability can occur in the southern region of HD 142527, whereas Am and β are close to the boundary between active and inactive parameters for the magnetorotational instability³¹. This is consistent with near-infrared observations, which suggest that the HD 142527 disk is turbulent³². Furthermore, the magnetic field strength of 0.3 mG at ~200 au from the central star is consistent with a prediction from numerical simulations³³ and within the upper limit constrained by Zeeman observations³⁴. By assuming that B_z has a power-law profile B_z ∝ r^p, the power-law index was estimated to be p ≈ −1 from the relative strengths of B_r/B_z (Methods). Such a power-law profile of the vertical magnetic field suggests that the magnetic field increases moderately towards the centre. That the power exponent of the magnetic field is flatter than 2 suggests that the frozen magnetic flux has thawed and that both magnetic field advection and outward diffusion play a role^35,36.

The relative strengths of the three components of the magnetic field also have important implications for disk evolution processes^25,26,37. For example, angular momentum transport due to the magnetic field can be inferred from the Maxwell’s stress M, as ({M}_{rphi }={{B}_{r}{B}_{phi }}/{4uppi }) and ({M}_{zphi }={{B}_{z}{B}_{phi }}/{4uppi }) along the radial (mass accretion) and vertical (disk wind) directions, respectively. The alpha value is commonly used to evaluate the angular momentum transport as ({alpha }_{rphi }={{M}_{rphi }}/{rho {c}_{{rm{s}}}^{2}}={2{B}_{r}}/{{beta }_{phi }{B}_{phi }}) and ({alpha }_{zphi }={{M}_{zphi }}/{rho {c}_{{rm{s}}}^{2}}={2{B}_{z}}/{{beta }_{phi }{B}_{phi }}). In our analysis, β was found to be ~2.0 × 10². Therefore, the alpha values are ({alpha }_{rphi } approx{0.53}/{beta }=2.5times 10^{-3}) and ({alpha }_{zphi } approx{0.6}/{beta }=2.2times 10^{-3}), respectively. We, therefore, conclude that the angular momentum transports in the radial and vertical directions are comparable in the ~200 au region of the protoplanetary disk around HD 142527.

In this paper, we demonstrated that dust polarization can be used as a magnetic tracer for protoplanetary disks when the dust grains are as small as those in the interstellar medium. The three-dimensional magnetic field can be probed for a disk with moderate inclination. Our analyses allowed us to estimate several physical parameters related to the magnetic field in the disk and provide insights into angular momentum transfer, which is an important process determining the evolution of the disk. This method is also applicable to other disks. Thus, our observations shed light on the study of the magnetic field in protoplanetary disks.

Observations and data reduction

We used archival data from ALMA projects 2015.1.00425.S, 2017.1.00987.S, 2018.1.01172.S and 2022.1.00406.S for this study. The results of the projects 2015.1.00425.S and 2018.1.01172.S have been presented previously^10,12,13. These projects observed HD 142527 with full Stokes polarization. The basic parameters of these observations are summarized in Supplementary Table 1. The polarization calibrator was observed 3 or 4 times with ~6 min of integration time during each observing schedule to calibrate the instrumental polarization (D-terms), cross-hand delay and cross-hand phase. The reduction and calibration of the data were performed in a standard manner using CASA³⁸.

Project 2015.1.00425.S was carried out on 11 March 2016 during its cycle 3 operation and on 21 May 2017 during its cycle 4 operation. The correlator set-up consisted of four spectral windows with a bandwidth of 1.75 GHz centred at sky frequencies of 336.5, 338.5, 348.5 and 350.5 GHz, providing a total bandwidth of ~7.5 GHz. The total integration times for the target were 73 min in cycle 3 operation and 80 min in cycle 4 operation.

Project 2017.1.00987.S was carried out on 7 January 2018 during its cycle 5 operation. The correlator set-up consisted of four 2 GHz spectral windows centred at sky frequencies of 137.995, 139.932, 149.995 and 151.995 GHz, providing a total bandwidth of ~8 GHz.

Project 2018.1.01172.S was carried out on 29 April 2019 during its cycle 6 operation. The correlator set-up consisted of four spectral windows. One band was set for 1.3 mm dust continuum emission centred at sky frequency 234.5 GHz with a bandwidth of 2 GHz. The other three windows had a bandwidth of 59 MHz and were centred on CO (J = 2 → 1), ¹³CO (J = 2 → 1) and C¹⁸O (J = 2 → 1) emission lines. In this paper, we used the spectral window of the 1.3 mm dust continuum emission.

Project 2022.1.00406.S was carried out on 26, 27, 30 and 31 March 2023 during its cycle 9 operation. The correlator set-up consisted of four spectral windows, and one band was set for 2.7 mm dust continuum emission centred at sky frequency 112 GHz with a bandwidth of 2 GHz. The other three windows had a bandwidth of 59 MHz and were centred on ¹³CO (J = 1 → 0), C¹⁸O (J = 1 → 0) and CN (N = 1 → 0, J = 1/2 → 1/2 and F = 1/2 → 1/2) emission lines. In this paper, we used the spectral window of the 2.7 mm dust continuum emission.

Imaging

Stokes I, Q and U images were generated from the calibrated visibility data by the CASA task tclean. To improve image fidelity, we performed an iterative phase-only self-calibration using the initial clean image as the first model image in CASA 6.5.1 for all projects except 2022.1.00406.S, as it was not improved by self-calibration. The time needed to solve the complex gain was reduced from infinity to 300 s and finally to 60 s. All images were generated with Briggs weighting. The robust parameters and the resulting parameters for the images are summarized in Supplementary Table 2. The PI was calculated from Stokes Q and U and has a positive bias because the Stokes Q and U components give the polarized intensity as (sqrt{{Q}^{2}+{U}^{,2}}). This bias is particularly noticeable in low signal-to-noise measurements. We, therefore, debiased the PI image as ({rm{PI}}=sqrt{{Q}^{2}+{U}^{2}-{sigma }_{{rm{PI}}}^{2}}), where σ_PI is the root-mean-square noise level^39,40. The PF image was generated by PF(%) = PI/(Stokes I) × 100 when the PI emission was above the noise level 3σ_PI. The polarization vectors were calculated from Stokes Q and U as ({rm{PA}}=frac{1}{2}arctanleft(frac{{rm{Stokes}},U}{{rm{Stokes}},Q}right)) when the PI emission was above the 3σ_PI noise level.

Imaging the differences between polarization vectors and polarization fractions

To avoid confounding by different spatial frequency components, the images of the polarization vector differences (Supplementary Fig. 1) were generated using the same uv distance. Note that uv is defined as the projected length of the baseline vector between two antennas in the Fourier (uv) plane. The uv distance is given by √u²+v², and it corresponds to the spatial frequency sampled by the interferometer. For ΔPA_{1.3mm−0.87mm} = PA_1.3mm − PA_0.87mm, the Stokes Q and U images for the 0.87 and 1.3 mm polarization data were generated using uv distances of 1.4 × 10⁴λ to 5.5 × 10⁵λ. The Stokes Q and U images for the 0.87 mm data were smoothed to a beam size of (0.51{7}^{{primeprime} }times 0.45{2}^{{primeprime} }) with a position angle of 64.2° to match the 1.3 mm data. The polarization vectors were then calculated from the Stokes Q and U images when the PI emission was above the 3σ levels. The σ_PI was derived to be 35 and 23 μJy per beam for 0.87 and 1.3 mm, respectively.

For ΔPA_{2.1mm−0.87mm} = PA_2.1mm − PA_0.87mm, the Stokes Q and U images for the 0.87 and 1.3 mm polarization data were generated using a uv distance of 1.4 × 10⁴λ to 1.2 × 10⁶λ. The Stokes Q and U images for the 0.87 mm data were smoothed to a beam size of (0.47{9}^{{primeprime} }times 0.40{5}^{{primeprime} }) with a position angle of 76.1° to match the 2.1 mm data. Then, the polarization vectors were calculated from the Stokes Q and U images. The polarization vectors were calculated when the PI emission was above the 3σ levels, where σ_PI was derived to be 33 and 7.6 μJy per beam at 0.87 and 2.1 mm, respectively. To make the images of the ratios of the polarization fractions (Supplementary Fig. 3), the polarization fractions were also calculated by the above method. The polarization fractions were calculated when the PI emission was above the 3σ levels.

For ΔPA_{2.7mm−0.87mm} = PA_2.7mm − PA_0.87mm, the Stokes Q and U images for the 0.87 and 2.7 mm polarization data were generated using a uv distance of 1.4 × 10⁴λ to 4.7 × 10⁵λ. The Stokes Q and U images for the 0.87 mm data were smoothed to a beam size of (0.71{1}^{{primeprime} }times 0.67{8}^{{primeprime} }) with a position angle of 50.6° to match the 2.7 mm data. Then, the polarization vectors were calculated from the Stokes Q and U images. The polarization vectors were calculated when the PI emission was above the 3σ levels, where σ_PI was derived to be 54 and 7.4 μJy per beam at 0.87 and 2.7 mm, respectively. To make the images of the ratios of the polarization fractions (Supplementary Fig. 3), the polarization fractions were also calculated with the above method. The polarization fractions were calculated when the PI emission was above the 3σ levels.

To evaluate the wavelength dependence of the polarization patterns, we calculated the differences for the polarization vectors at 1.3, 2.1 and 2.7 mm with respect to the 0.87 mm polarization vectors by deriving ΔPA_{1.3mm−0.87mm} = PA_1.3mm − PA_0.87mm, ΔPA_{2.1mm−0.87mm} = PA_2.1mm − PA_0.87mm and ΔPA_{2.7mm−0.87mm} = PA_2.7mm − PA_0.87mm. Supplementary Fig. 1 shows the differences in the polarization vectors. In the northern region, ΔPA_{1.3mm−0.87mm}, ΔPA_{2.1mm−0.87mm} and ΔPA_{2.7mm−0.87mm} are all ~90° on the ridge of the horseshoe structure, indicating that PA_0.87mm is perpendicular to PA_1.3mm, PA_2.1mm and PA_2.7mm. By contrast, in the southern region, the polarization vectors are almost in the same direction (ΔPA_{1.3mm−0.87mm}, ΔPA_{2.1mm−0.87mm} and ΔPA_{2.7mm−0.87mm} ≈ 0°). Supplementary Fig. 2 plots ΔPA_{1.3mm−0.87mm}, ΔPA_{2.1mm−0.87mm} and ΔPA_{2.7mm−0.87mm} as functions of position angles of 110–250° on the ridge of the southern part of the horseshoe structure. The plotted data were taken from the pixels at the ridge position with a Nyquist sampling of 11° for ΔPA_{1.3mm−0.87mm}, 10° for ΔPA_{2.1mm−0.87mm} and 15° for ΔPA_{2.7mm−0.87mm}. The differences for the polarization vectors were derived to be ΔPA_{1.3mm−0.87mm} = 3.9 ± 0.6°, ΔPA_{2.1mm−0.87mm} = 1 ± 1° and ΔPA_{2.7mm−0.87mm} = −5.2 ± 3.4°. Supplementary Figs. 1 and 2 show that the polarization vectors have almost the same directions across all the wavelengths. Although larger angular offsets are seen at longer wavelength separations, this may be due to the weaker emission at longer wavelengths. The 2.7 mm polarization emission in band 3 is the weakest emission with the largest beam size, leading to larger uncertainties in the polarization vectors due to the noise. In addition to the polarization vectors, Supplementary Fig. 3 shows the ratios of the polarization fraction P_frac,1.3mm/P_frac,0.87mm, P_frac,2.1mm/P_frac,0.87mm and P_frac,2.7mm/P_frac,0.87mm (Methods). In the northern region, the polarization fractions differ significantly, with P_frac,1.3mm/P_frac,0.87mm reaching values as high as ~5 and P_frac,2.1mm/P_frac,0.87mm and P_frac,2.7mm/P_frac,0.87mm reaching values as high as ~10. By contrast, in the southern region, the polarization fractions are almost the same values, as all the polarization fraction ratios are ~1.

The dual-layer model for the angular offset

In Fig. 3, we have considered the single-layer magnetic field to estimate the relative strengths. However, B_ϕ and B_r were assumed to be in opposite directions in the upper and lower disk layers because these magnetic fields are generated by the entrainment and accretion of the initial B_z due to the rotation (Extended Data Fig. 1). Thus, the observed offset may be due to a combination of the upper disk layer (Δϕ₊) and the lower disk layer (Δϕ₋).

By assuming that the radial and toroidal magnetic fields in the lower layer have the same strengths but in opposite directions to those in the upper layer, the angular offset of the lower layer (Δϕ₋) was estimated as (Delta {phi }_{-}=arctanscriptstylefrac{-{({B}_{r})}_{phi perp }+{({B}_{z,mathrm{proj}})}_{phi perp }}{-{B}_{phi }-{({B}_{r})}_{phi }+{({B}_{z,mathrm{proj}})}_{phi }}). Because the 0.87 mm dust continuum emission is optically thin (τ < 0.1), the angular offset (Δϕ) should be the weighted average of Δϕ₊ and Δϕ₋:

where the subscripts Δϕ₊ and Δϕ₋ in the Stokes Q and Stokes U parameters represent emission from the upper and lower layers, respectively. The (overline{Delta phi }) value is not simply determined by the combination of Δϕ₊ and Δϕ₋ but by the combination of the Stokes Q and U emission.

The polarization fraction was calculated with (papprox {p}_{mathrm{max}}{mathrm{sin}}^{2}theta), where θ is the projected length of the elongated dust grains along the plane of the sky^41,42 and ({p}_{max }) is defined as the intrinsic polarization, which is the maximum polarization when the grains are perfectly aligned in the same direction. However, ({p}_{max }) can be ignored in this calculation because we assume the same ({p}_{max }) between the upper and lower layers. Then, the (overline{Delta} {phi }) value is independent of ({p}_{max }). Therefore, Stokes Q and Stokes U were derived only from ({rm{sin}}^{2}theta) in both layers.

Estimating the magnetic field strength and power-law of the vertical magnetic field

From the relative strengths of the three-dimensional magnetic fields, we estimated the magnetic field strength and the power-law profile of the vertical magnetic field with radius. By assuming a balance between the generation of B_ϕ by vertical shear motion and the dissipation of B_ϕ by ambipolar diffusion, B_ϕ can be estimated⁴³ as

where H, r, η_A, v_ϕ and Ω_K are the gas scale height, radius, ambipolar resistivity, Kepler velocity and Kepler frequency. Therefore, the ambipolar resistivity can be written as ({eta }_{{rm{A}}}=scriptstylefrac{{H}^{3}{varOmega }_{{rm{K}}}}{r}frac{{B}_{z}}{{B}_{phi }}). The gas scale height is given by H = c_s/Ω_K, where c_s is the speed of sound. Here we assumed that the dust scale height is the same as the gas scale height because it has been suggested that this disk is highly turbulent³². The turbulent strength (α) and dust size (a_dust) were estimated to be α ≈ 0.3 and a_dust ≈ 3 μm, respectively, from near-infrared observations with the compact dust assumption. Then, the dust scale height can be calculated as ({H}_{{rm{dust}}}={left(1+frac{{rm{St}}}{alpha }frac{1+2{rm{St}}}{1+{rm{St}}}right)}^{-1/2}{H}_{{rm{gas}}}) (refs. ^44,45), where St = (πρa_dust)/(2Σ_gas) is the Stokes number estimated from the dust internal density (ρ), dust size (a_dust), and gas surface density (Σ_gas). Therefore, H_dust ≈ H_gas when the dust size is less than 100 μm, assuming a gas surface density of 0.2 g cm⁻² and a turbulent strength α = 0.3. Note that the turbulent strength has been studied for different porous dust grains and dust surface densities by keeping a fixed value of St/α (ref. ³²). Even if the dust structure and dust surface density are changed, our assumption of H_dust ≈ H_gas is robust as long as St/α does not change. In addition, we assumed that the relative strengths of the magnetic field do not vary much within the dust scale height. This is suggested to be a reasonable assumption for B_z, because the variation of B_z within the gas scale height is estimated to be (Delta {B}_{z}approx({H}/{r}){B}_{r}approx({H}{r})/{B}_{z}) (ref. ⁴⁶). Here, r ≈ 200 au and H ≈ 20 au. Therefore, ΔB_z can be negligible. By contrast, the variations of B_r and B_ϕ within the gas scale height are unknown parameters and have been studied by magnetohydrodynamics simulations. If the single-layer model applies, it can be speculated that the strengths of B_r and B_ϕ do not vary as much within the scale height²⁹. In the dual-layer model, the polarized emission is from the relatively upper layer because B_ϕ will fall at the disk midplane so that its direction can reverse, whereas the observations indicate that B_ϕ is the dominant component. Note that the emission was estimated to be optically thin (τ ≈ 0.1) at a wavelength of 0.87 mm, suggesting that the polarized emission is the average value along the vertical direction of the scale height and could mainly trace the midplane region. This also indicates that the single-layer model is better than the dual-layer model.

By considering a central star of mass 2.4 M_⊙, temperature 30 K and ({{B}_{z}}/{{B}_{phi }}) = 0.23, η_A was derived to be 2.3 × 10¹⁷ cm² s⁻¹. Then, because the gas density was small enough, the ambipolar resistivity was approximated as η_A = B²/(4πCγρ^3/2) (refs. ^43,47). Here C is (C=scriptstylesqrt{{{m}_mathrm{i}^{2}{zeta }_{{rm{CR}}}}/{{m}_mathrm{g}{beta }_mathrm{r}}}), where m_i and m_g are the masses of ions and neutral particles. We adopted m_i = 29m_H and m_g = 2.34m_H assuming that the major ion is HCO⁺, where m_H is the hydrogen mass. The HCO⁺ molecule has been suggested to be the most abundant ion when the abundance ratio of n(CO)/n(electron) is higher than ~10³ (ref. ⁴⁸). Although the electron abundance is not well constrained, it is usually assumed to be X(e) ≲ 10⁻⁹. On the other hand, the CO abundance has been estimated to be X(CO) ≈ 10⁻⁴ for the HD 142527 disk¹¹. Therefore, we assumed that the HCO⁺ molecule is the most abundant ion for this region. ζ_CR is the cosmic ray ionization rate, which we assumed was 10⁻¹⁷ s⁻¹. β_r is the recombination rate and was assumed to be ({beta }_mathrm{r}={beta }_mathrm{r,0}{({T}/{300,{rm{K}}})}^{-0.67}), where β_r,0 = 2.4 × 10⁻⁷ cm³ s⁻¹. (gamma ={{langle sigma vrangle }_{{rm{in}}}}/{({m}_mathrm{g}+{m}_mathrm{i})}), where 〈σv〉_in is the rate coefficient for collisional momentum transfer between ions and neutrals. We assumed 〈σv〉_in = 1.3 × 10⁻⁹ cm³ s⁻¹. We estimated the magnetic field strength with this formula. The gas density ρ was derived to be ρ = Σ_g/H, where Σ_g is the gas surface density, which was derived to be Σ_g = 0.2 g cm⁻² (ref. ¹¹). From these parameters, the magnetic field strength was finally derived to be B ≈ 0.3 mG. The Am parameter was also estimated from the magnetic field strength. The Am parameter is defined as ({rm{Am}}={{v}_{{{rm{A}}}^{2}}}/{{eta }_{{rm{A}}}{varOmega }_{{rm{K}}}}={{B}^{2}}/({{eta }_{{rm{A}}}{varOmega }_{{rm{K}}}times 4uppi rho})), where v_A is the Alfvén velocity, ({v}_{{rm{A}}}=sqrt{{{B}^{2}}/{4uppi rho }}). By assuming B ≈ 0.3 mG, η_A = 2.3 × 10¹⁷ cm² s⁻¹ and ρ = 6.6 × 10⁻¹⁶ g cm⁻³, the Am parameter was derived to be Am ≈ 0.4. The plasma beta is defined as (beta ={8uppi rho {c}_{{rm{s}}}^{2}}/{{B}^{2}}) and was derived to be β ≈ 2.0 × 10².

By assuming that B_z has a power-law profile B_z ∝ r^p, the power-law index p was estimated from the relative strengths between B_r and B_z. The radial magnetic field B_r was obtained approximately as ({B}_{r}={varPsi (r)}/{2uppi {r}^{2}}), where Ψ(r) is the magnetic flux enclosed within a radius r and is given by (varPsi =int_{0}^{r}2uppi r{B}_{z},mathrm{d}r). By assuming the power-law profile ({B}_{z}={B}_{z,0}{({r}/{1,{rm{au}}})}^{p}), (Psi (r)=scriptstyleint_{0}^{r}2uppi r{B}_{z,0}{left(frac{r}{1,{rm{au}}}right)}^{p},mathrm{d}r=frac{2uppi {B}_{z,0}}{p+2}{left(frac{r}{1,{rm{au}}}right)}^{p+2}), where B_z,0 is the vertical magnetic field strength at a radius of 1 au. From the magnetic flux, B_r can be written as ({B}_{r}=scriptstylefrac{2uppi {B}_{z,0}}{p+2}{left(frac{r}{1,{rm{au}}}right)}^{(p+2)}frac{1}{2uppi {r}^{2}}=frac{{B}_{z,0}}{(p+2)}{left(frac{r}{1,{rm{au}}}right)}^{p}). Therefore, (frac{{B}_{r}}{{B}_{z}}=frac{1}{(p+2)}). As (frac{| {B}_{r}| }{| {B}_{z}| }approx 1) around a radius of 200 au, the power-law index was derived to be p ≈ −1. This implies that the magnetic field increases moderately towards the centre. That the power exponent of the magnetic field is flatter than 2 suggests that the frozen magnetic flux has thawed and that both magnetic field advection and outward diffusion both play a role^35,36.