The C IV doublet in the UV has long been associated with accretion in T Tauri stars. However it is still unclear where and how the lines are formed. Here we present a new C IV line model based on the currently available accretion shock and accretion flow models. We assume axisymmetric dipolar accretion flows with different energy fluxes and calculate the properties of the accretion shock. We use Cloudy to obtain the carbon level populations and calculate the emerging line profiles assuming a planeparallel geometry near the shock. Our model generally reproduces the intensities and shapes of the C IV emission lines observed from T Tauri stars. We find that the narrow component is optically thin and originates in the postshock while the broad component is optically thick and emerges from the preshock. We apply our model to seven T Tauri stars from the Hubble Ultraviolet Legacy Library of Young Stars as Essential Standards Directors Discretionary program ULLYSES for which the consistently determined accretion shock properties are available. We can reproduce the observations of four stars finding that the accretion flows are carbondepleted. We also find that the chromospheric emission accounts for less than 10 of the observed C IV line flux in accreting T Tauri stars. This work paves the way toward a better understanding of hot line formation and provides a potential probe of abundances in the inner disk.
['T Tauri stars', 'Stellar accretion', 'Ultraviolet astronomy', 'Protoplanetary disks', 'Astronomical models', 'Radiative transfer', 'Chemical abundances', '1681', '1578', '1736', '1300', '86', '1335', '224', 'Astrophysics - Solar and Stellar Astrophysics', 'Astrophysics - Earth and Planetary Astrophysics']
A Model of the C IV  1548 1550 Doublet Line in T Tauri Stars
{'A Model of the C IV λλ 1548, 1550 Doublet Line in T Tauri Stars': '<!-- image -->  \n1 Institute for Astrophysical Research, Department of Astronomy, Boston University, 725 Commonwealth Ave., Boston, MA 02215, USA 2 Department of Physics & Astronomy, Amherst College, C025 Science Center 25 East Drive, Amherst, MA 01002, USA  \n3 Department of Astronomy, University of Michigan, 1085 South University Ave., Ann Arbor, MI 48109, USA  \n4 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  \n5 Laboratory for Atmospheric and Space Physics, University of Colorado Boulder, Boulder, CO 80303, USA  \n6 Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, 77 Massachusetts Ave., Cambridge, MA 02139, USA  \n7 Max-Planck-Institut fur Astrophysik, Karl-Schwarzschild-Stra β e 1, 85748 Garching b. Munchen, Germany 8 Hamburger Sternwarte, Gojenbergsweg 112, D-21029, Hamburg, Germany  \n(Received 2024 June 21; Revised 2024 September 13; Accepted 2024 September 13)  \nSubmitted to ApJ', 'ABSTRACT': "The C iv doublet in the UV has long been associated with accretion in T Tauri stars. However, it is still unclear where and how the lines are formed. Here, we present a new C iv line model based on the currently available accretion shock and accretion flow models. We assume axisymmetric, dipolar accretion flows with different energy fluxes and calculate the properties of the accretion shock. We use Cloudy to obtain the carbon level populations and calculate the emerging line profiles assuming a plane-parallel geometry near the shock. Our model generally reproduces the intensities and shapes of the C iv emission lines observed from T Tauri stars. We find that the narrow component is optically thin and originates in the postshock, while the broad component is optically thick and emerges from the preshock. We apply our model to seven T Tauri stars from the Hubble Ultraviolet Legacy Library of Young Stars as Essential Standards Director's Discretionary program (ULLYSES), for which consistently determined accretion shock properties are available. We can reproduce the observations of four stars, finding that the accretion flows are carbon-depleted. We also find that the chromospheric emission accounts for less than 10 percent of the observed C iv line flux in accreting T Tauri stars. This work paves the way toward a better understanding of hot line formation and provides a potential probe of abundances in the inner disk.  \nKeywords: T Tauri stars (1681), Stellar accretion (1578), Ultraviolet astronomy (1736), Protoplanetary disks (1300), Astronomical models (86), Radiative transfer (1335), Chemical abundances (224)", '1. INTRODUCTION': "Magnetospheric accretion has been the paradigm of accretion in low-mass, pre-main sequence stars for more than two decades, with support from theoretical and observational studies. Numerical models of magnetospheric accretion in young pre-main sequence stars (T  \nCorresponding author: Thanawuth Thanathibodee  \nthanathi@bu.edu  \nTauri stars; TTS) have been developed to explain observations in hydrogen lines to great length (Hartmann et al. 1994; Muzerolle et al. 1998, 2001; Bouvier et al. 2007; Kurosawa & Romanova 2013; Esau et al. 2014; Tessore et al. 2021). Applying these models has provided insight into the overall structure of the magnetospheric flow, barring some assumptions. Similarly, advances in infrared interferometric observations of the closest T Tauri stars (e.g., Gravity Collaboration et al. 2020, 2023a,b) have confirmed this picture of accretion onto young stars.  \nAn accretion shock on the stellar surface is one of the consequences of magnetospheric accretion in T Tauri stars (see Hartmann et al. 2016, for review). A shock is formed as the accretion flow material reaches the photosphere near free-fall velocities (hundreds of km s -1 ). The kinetic energy from the flow is reprocessed and heats the region before the shock (the preshock), immediately after the shock (the postshock), the underlying photosphere (the heated photosphere), as well as its immediate surroundings (Calvet & Gullbring 1998; Colombo et al. 2019). Models of the accretion shock have been used to explain the excess continuum emission in the near-UV of accreting TTS (Calvet & Gullbring 1998). Variations of such models have been used to infer the structure of the accretion shock hotspots (Ingleby et al. 2013; Espaillat et al. 2021; Robinson et al. 2022; Pittman et al. 2022), suggesting the density stratification at the base of the accretion flow. Such stratification is also inferred from X-ray observations (e.g., Schneider et al. 2018).  \nIn addition to continuum excess emission, the regions associated with the accretion shock are also thought to produce emission lines. However, where and how observed lines form is still uncertain. Much progress has been made in using X-ray emission lines to improve our understanding of properties of accretion shock regions, such as densities, temperature, geometry, and chemical abundances (e.g., Kastner et al. 2002; Stelzer & Schmitt 2004; Drake et al. 2005; Argiroffi et al. 2005; Brickhouse et al. 2010; Argiroffi et al. 2017). Shockoriginated emission lines are also found in the optical, such as He I λ 5876, which is used to determine magnetic obliquity in TTS (McGinnis et al. 2020).  \nAtomic emission lines originating in the magnetospheric flow from the inner disk onto the star also provide an indirect probe of the inner disk composition (Micolta et al. 2023). Carbon is an especially critical tracer of this region, as it is one of the most abundant elements in the Universe after hydrogen and helium. It is a key element of (Earth-like) life, and carbon-bearing species play important roles in the chemistry of protoplanetary disks (e.g., Oberg et al. 2023). Therefore, detecting and characterizing carbon-bearing species can provide many insights into the properties of protoplanetary disks, including their chemical composition, density, and temperature structures. Carbon-bearing molecules have been studied extensively in the sub-mm and infrared using facilities such as ALMA, the VLT, and JWST. These lines form in the colder regions of protoplanetary disks. Probing the inner disk is more challenging. The C I atomic lines in the near-infrared (e.g., 0.908, 0.94, 1.069, 1.176, 1.45, 2.166 µ m) provide a way, but they have been observed in only a handful of sources (McClure 2019; McClure et al. 2020).  \nAnother way to indirectly probe the inner disk is by studying the magnetospheric flow from the disk onto the star. Since the flow's temperatures are high, all  \nelements are expected to be atomic and possibly ionized. Therefore, studying atomic emission lines formed in the accretion flow or the shock can be used to determine the elemental compositions of the inner disk (Micolta et al. 2023).  \nThe C IV lines at 1548 and 1550 ˚ A are some of the most prominent lines in the spectra of T Tauri stars (see review by Schneider et al. 2020), and they were recognized very early to be related to accretion (Ardila et al. 2002; Johns-Krull et al. 2000; Valenti et al. 2000; Calvet et al. 1996; Valenti et al. 1993). Each of the lines in the doublet is generally composed of a narrow component (NC) with low redshifted velocity (a few 10 km s -1 ) and a broad component (BC) with velocity widths on the order of a few hundred km s -1 and a wide range of the centroid velocity ( ∼ -20 to ∼ +200 kms -1 ; Ardila et al. 2013). Using extensive data from the HST DAO program (PI: G. Herczeg, PID: 11616), Ardila et al. (2013) analyzed the C IV lines of 28 TTSs and confirmed previous studies that the line luminosity correlates with mass accretion rates, suggesting that the line forms in the accretion process. The line profiles, on the other hand, are harder to explain. The Gaussian decomposition results of Ardila et al. (2013) suggest that only half of the CTTS in their sample are consistent with the scenario where the BC is formed in the preshock, as the centroid velocity of the BC is less than four times that of the NC, as predicted by the jump condition for the shock. However, as we will discuss later, we note that their interpretation does not consider geometry effects, such that the emission from the preshock may peak at lower velocities.  \nSeveral theoretical works have attempted to explain the observed line profiles to study how the line forms. Lamzin (2003a,b) calculated the properties of the accretion shock and the C IV lines using axially symmetric shock and radial accretion and found that the model profiles are inconsistent with the observations. Colombo et al. (2016) synthesized the C IV line from MHD simulations with fragmented accretion streams and trains of accretion blobs and found a reasonable agreement between the synthetic line profiles and the observed profiles, but the line flux levels are one to two orders of magnitude lower than that of the observations. In this work, we expand upon the work of Lamzin (2003b) and present a new model of the C IV line formation based on the accretion shock model of Calvet & Gullbring (1998) and the accretion flow model of Muzerolle et al. (2001). The goals of this paper are twofold: (1) to understand where and how the C IV line forms in the context of the accretion shock and (2) to use the C IV line to infer properties of the innermost disk regions around accreting T Tauri stars. We describe our model in Section 2 and provide some model results in Section 3. In Section 4, we apply the model to a subset of ULLYSES targets  \nin Orion OB1 and discuss the results in Section 5. We provide conclusions and a summary in Section 6.", '2. DESCRIPTION OF THE MODEL': 'As in previous studies (Ardila et al. 2013; Lamzin 2003a,b), we consider that the C IV doublet forms near the shock created by accreting material reaching the photosphere. In this section, we provide a detailed description of the model and the line calculation.', '2.1. Geometry': "The large-scale geometry of the model is illustrated in Figure 1. A given point on the star's surface is defined using spherical coordinates ( θ , ϕ ) centered at the star. θ is measured from the pole, and ϕ is defined such that the vector pointing toward the Earth is ( θ , ϕ )=( i , π/ 2), where i is the star's inclination. We assume that the shock is created at the base of the axisymmetric accretion flow in a dipole geometry, which lands perpendicular to the surface. As a result, the shock surface is an annulus on the star. Within this framework, a filling factor f ' subtends an area on the stellar surface between any two co-latitudes θ 1 and θ 2 as  \nf ' = 1 2 [cos( θ 1 ) -cos( θ 2 )] . (1)  \nSimilarly to Pittman et al. (2022), we assume that the ring is composed of multiple flows at three different energy fluxes: high energy flux F h , medium energy flux F m , and low energy flux F l . Motivated by results from numerical simulations and observations (e.g., Romanova et al. 2004; Espaillat et al. 2021), we assume that two medium-energy flux flows blanket the central flow, which has the highest energy flux (and highest density), and together they are sandwiched by two low-energy flux flows. Therefore, there are five rings in total for each hemisphere (c.f., Figure 1). The same geometry also occurs in the other hemisphere since we expect that the material in the accretion flows lands on both hemispheres. That is  \nf h,m,l = 2 f ' h,m,l , (2)  \nwhere f is the total filling factor of a given energy flux. The co-latitude for each ring is defined by two associated magnetic radii R in and R out through the dipolar flow equation as  \nθ = 1 2 ( θ in + θ out ) (3)  \n= 1 2 [ sin -1 ( √ R ∗ R in ) +sin -1 ( √ R ∗ R out )] . (4)  \nFigure 1 shows the large-scale geometry of the multicomponent flow. We assume that the flows at different energy fluxes are independent. We also assume that the surrounding, large-scale accretion flows ( ∼ few stellar  \nradii) do not absorb any emission emerging from the shock.  \nAs in Calvet & Gullbring (1998), the shock at the star's surface is composed of the preshock and the postshock region. For simplicity, we assume that the flow lands perpendicular to the star's surface and that the shock structure's height ( ∼ the shock cooling distance) is much smaller than the stellar radius. Therefore, we assume plane-parallel geometry for the following radiative transfer calculation.  \nFor a given spot ( θ, ϕ ) on the surface, it can be shown (e.g., Gunther & Wawrzyn 2011) that the angle between the line of sight and the flow is given by  \ncos θ 0 = sin θ sin ϕ sin( i ) + cos θ cos( i ) , (5)  \nfor the upper hemisphere ( θ < π/ 2), and  \ncos θ 0 = sin θ sin ϕ sin( i ) -cos θ cos( i ) , (6)  \nfor the lower hemisphere ( θ > π/ 2). We only include in our calculation the point ( θ , ϕ ) where θ 0 < π/ 2, since the flow would be behind the star for a larger angle. This requirement limits the longitude for the ring to be seen between ϕ w and ϕ e .", '2.2.1. Initial Calculations': 'We use the 1D model of Calvet & Gullbring (1998), as updated by Robinson & Espaillat (2019), to calculate the structure of the preshock and postshock, and we refer the reader to those papers for details. The postshock density and temperature are calculated by solving the fluid equations coupled with hydrostatic equilibrium equations using emissivities from Cloudy (Ferland et al. 2017). The resulting emergent spectrum from the postshock F post is then used as an input to calculate the emergent spectrum of the preshock F pre , using Cloudy. These spectra are used, without modification, for subsequent runs to calculate ionization, continuum opacities, and carbon level populations. Figure 2 shows the geometry of the shock and the level population calculation.', '2.2.2. Subsequent Calculations': "After the initial run of the model to generate F post and F pre , we use Cloudy to calculate the relevant properties of the post and preshock for line radiative transfer. Cloudy is generally used to determine the steadystate properties of a slab of material irradiated by an external source given a density structure and, optionally, a temperature structure. In our model, the postshock is irradiated by the emission from the preshock ( F pre ); similarly, the preshock is irradiated by the postshock ( F post ). Some previously calculated properties of the slabs are fixed. The density and temperature structure of the postshock are fixed from the shock model calculation. For the preshock, we assume a constant density; Cloudy calculates the temperature along  \nFigure 1. Left: The overall geometry of the model. The star's equatorial plane is the same as that of the disk. Axisymmetric accretion funnels originate from the disk and land perpendicular to the star's surface. The vector toward the observer is defined to be in the yz plane with the colatitude i , i.e., the inclination of the star. Right: A zoom-in, side-view illustration of the large-scale geometry of accretion flows landing on the surface of a star. The darker color corresponds to the highest energy flux ( F h ). The post/preshock are small regions at the base of the flows. Gray dotted lines are the center of the flow used to define the colatitude θ for a given flow.  \n<!-- image -->  \nFigure 2. The schematic view of an accretion column and its associated geometry. An observer is viewing the flow at an angle θ 0 .  \n<!-- image -->  \nwith other properties. In both the postshock and the preshock, we also include turbulent velocities v t,pre/post of the flow and carbon abundance in the solar scale [ C/H ] = log( N C /N H ) ⋆ -log( N C /N H ) ⊙ , as parameters. In these Cloudy runs, we record the level population and other properties as a function of depth from the shock.", '2.3. Line Profile Calculation': "Since the post and preshock heights are smaller than the stellar radius, we assume plane-parallel geometry to calculate the line profiles. The overall flux is given by,  \nF ν = R 2 ∗ d 2 ∫ ϕ W ϕ E ( I ν,shock + I ν,⋆ ) cos θ 0 sin θ ∆ θdϕ, (7)  \nwhere we assume that the ring is thin and ∆ θ = ( θ in -θ out ), R ∗ is the radius of the star, and d is its distance. ϕ E and ϕ W refer to the easternmost and westernmost edges of the ring, at which point θ 0 = π/ 2. We assume a blackbody at the star's effective temperature for I ν,star . The θ 0 -dependent I ν,shock is composed of the  \ncontributions from the preshock and the postshock. We assume that the postshock is always viewed through the preshock, therefore  \nI ν,shock = I ν,pre + I ν,post e -τ ν,tot,pre , (8)  \nwhere τ ν,tot,pre is the total optical depth of the preshock. We follow the same steps to calculate the specific intensity and the optical depth of the preshock and the postshock, as shown below.  \nFor a given flow, the specific intensity at a position ( θ , ϕ ) depends on the angle between the accretion flow direction and the line of sight θ 0 (c.f., Figure 1). The plane-parallel assumption simplifies its calculation to the integration of the emissivity along the flow toward the star along r with a factor µ as  \nI ν ( µ ) = ∫ ( η l ( r ) + η c ( r )) e -µτ ν ( r ) µdr,  \nµ ≡ 1 cos θ ,  \n∞ r 0 (9) 0 (10)  \nwhere η l and η c are the line and continuum emissivity, respectively, and τ is the optical depth.  \nThe continuum emissivity is taken directly from Cloudy simulations. The line emissivity is  \nη ν,l = χ ν,l S l , (11)  \nwhere the line opacity is given by  \nχ ν,l = χ 0 ,l · V ν ( θ 0 , v t , v ( r )) (12)  \n= hν 0 4 π B lu N l ( 1 -N u g l N l g u ) · V ν ( θ 0 , v t , v ( r )) , (13)  \nwhere  \nB lu = g u g l B ul = g u g l c 2 2 hν 3 0 A ul , (14)  \nand line source function is given by  \nS l = 2 hν 3 0 c 2 [ N l g u N u g l -1 ] -1 . (15)  \nThe Voigt profile V ν depends on the turbulent velocity v t , assumed to be a constant parameter for a given flow, and the local properties, including the velocity v ( r ). It is calculated as  \nV ν = 1 π ∆ ν D H ( Γ r 4 π ∆ ν D , ν -ν ' ∆ ν D ) , (16)  \nwhere H is the Hjerting function and Γ r is the line radiative broadening parameter. Stark broadening is small for the electron densities considered here ( ∼ 10 12 -10 15 ; Elabidi et al. 2011), while the van der Waals broadening is expected to be minimal in the ionized environment that the line forms. The Doppler width is calculated from the thermal and turbulent velocities, ∆ ν D = ν 0 c √ v 2 term + v 2 t , and the shifted line-center frequency is ν ' = ν 0 (1 -v ( r ) cos θ 0 /c ). In the preshock, v ( r ) is a constant, adopted as a free-fall velocity from R in ; v pre = √ 2 GM ∗ (1 /R ∗ -1 /R in ) (Calvet & Gullbring 1998). In the postshock, v ( r ) decreases from 1 4 v pre ( r ) at the shock to ∼ 0 at the photosphere. v post ( r ) is obtained along with the postshock structure in Section 2.2.1.  \nThe optical depth along the flow is given by  \nτ ν ( r ) = ∫ ∞ r ( χ ν,l ( r ) + χ c ( r )) dr, (17)  \nwhere the continuum opacity χ c is computed by Cloudy.  \nAll relevant atomic parameters are retrieved from NIST (Kramida et al. 2023) and are summarized in Table 1.", '2.4. Model Limitations': "We make several approximations and assumptions to our model, and here we list these limitations. These assumptions are made while keeping the following goal in mind: to study how the geometry of the emitting region and the large-scale accretion properties affect the C IV line.  \nFirst, several works have shown that the accretion shock is dynamic, and time-dependent models may be relevant (e.g., de S'a et al. 2019). We note, however, that the timescale of quasi-periodic oscillation (QPO) cycles in the time-dependent models are on the order of 100-600 s. In contrast, the observational timescales (e.g., telescope exposure times) are on the order of 1000-3000 s or more. Dedicated analysis of TW Hya's observations with sufficient time resolution also found no evidence of this variability (Drake et al. 2009; Gunther et al. 2010, but, see also Siwak et al. 2018). Other MHD simulations of accretion shocks, such as those presented in (Orlando et al. 2010), show that these oscillations disappear in  \nTable 1. Atomic Parameters of the C IV lines.  \nthe 2D analysis. Therefore, a time-average model will be needed to explain the observation, and we opt to use steady-state approximation following most of the literature on this topic.  \nSecond, we used a readily available and widely utilized code, Cloudy, to calculate the population. This code uses the escape probability method for radiative transfer, and the results are not exact (Hubeny 2001). According to the Cloudy documentation Hazy 2 , this approximation may be problematic for transitions where the lower and upper levels are in the excited states, but it notes that resonance lines (e.g., C IV , with the lower level at electronic ground state) should be fine. Therefore, as a first approximation, we adopted Cloudy's formalism.  \nThird, we adopted a strictly 1D, plane-parallel geometry in calculating the flow structure and level population. That is, we do not consider the angular dependence of the irradiation in the steps described in Section 2.2, following Calvet & Gullbring (1998) and Robinson & Espaillat (2019). This omission may introduce some uncertainty in the level population and, subsequently, the line flux, especially near the edge of the flow. However, as pointed out by observations and simulations (e.g., Espaillat et al. 2021), the outer parts of the flow may have a lower density than the inner part, and this part contributes less to the overall flux/profile, as will be shown later in this work. Therefore, the 1D approximation is justified.", '3. MODEL RESULTS': 'In this section, we showcase the results using the models calculated for CVSO 104 (M ⋆ = 0 . 37 M ⊙ , R ⋆ = 1 . 66 R ⊙ ; Manara et al. 2021) and explore various properties of the flows. We also study the effects of varying turbulent velocities ( v t = 10 to 100 km s -1 ), carbon  \nFigure 3. The macroscopic and microscopic properties of an accretion flow with F = 10 11 erg s -1 cm -2 , M ⋆ = 0 . 37 M ⊙ , and R ⋆ = 1 . 66 R ⊙ . The horizontal axis indicates the distance from the transition between the heated photosphere and the postshock. The vertical dashed line marks the location of the shock. In the fourth row, the dashed lines indicate the fraction at each ionization state, assuming LTE, for comparison.  \n<!-- image -->  \nabundances ([C / H] = 0 to -3), inclinations ( i = 0 · to 90 · ), and truncation radii (3 . 14 and 5 R ⋆ ) on the C IV line properties and profiles.', '3.1. Properties of the flows': "After solving the density and velocity structure of the postshock and preshock and the temperature of the postshock (Section 2.2), we calculate the atomic properties related to the line profile calculation (see Section 2.2.) Figure 3 shows the structure of the flow and the resulting ionization and level population. As in Figure 2, the observer is viewing the flow from the right through the preshock, toward the shock, and finally, the postshock.  \nBy construction, the hydrogen number density N H and the flow velocity v ( r ) of the preshock are constant as the depth of the preshock is much smaller than the star's radius. After the shock, the temperature and density of the gas sharply increase while the speed suddenly drops, as expected. The fourth row of Figure 3 shows the ionization fraction of the carbon atom. In the post-  \nshock, C IV is only present in a thin, ∼ 10 3 cm layer near the transition to the heated photosphere. In contrast, C IV is the dominant ionization level for a much thicker ( ∼ 10 7 cm) region in the preshock.  \nThe C IV doublet lines are formed between the transitions of the first three electronic levels of C 3+ : [1] 2 S 1 / 2 , [2] 2 P 1 / 2 , and [3] 2 P 3 / 2 (c.f., Table 1). The last row of Figure 3 shows the percentage of each electronic level of the C IV ionization level. In every region, most electrons in the C IV atoms are in the ground state, as the other levels account for ∼ 1% of the population. In the top row of Figure 4, we compare the number density of the C IV electronic states across the shock region. While the density of the preshock is lower due to a higher velocity and constant mass flux, the number densities of C IV are comparable to the peak number density in the postshock and at a much larger extent. For comparison, we plot the LTE level population assuming the density and temperature of the given region, shown as dashed lines at the top of Figure 4. The number densities in the postshock depart slightly from LTE, with the densities  \nFigure 4. Atomic and radiative properties of the postshock and preshock flows with F = 10 11 erg s -1 cm -2 and M ⋆ = 0 . 37 M ⊙ , R ⋆ = 1 . 66 R ⊙ . The shock is marked with the vertical dashed line, and the distance is measured outwardly from the postshockheated photosphere transition. The observer is viewing through the flow from the right toward the star on the left. In the postshock, the dashed lines show the quantities with the LTE level populations. The LTE populations in the preshock are much lower than the range presented here and are not shown. In the second row, the dotted grey line is the Planck function, plotted here for comparison. The Planck function for the postshock traces the same curve as the LTE source function (dashed line). The units for S l , χ 0 ,l , and η 0 ,l are erg s -1 cm -2 sr -1 Hz -1 , cm -1 , and erg s -1 cm -3 sr -1 Hz -1 , respectively.  \n<!-- image -->  \npeaking at the same location at ∼ 10 2 . 5 cm. Closer to the shock, the C IV number densities decrease quickly as the atoms are ionized due to the high temperatures. The LTE population is much lower in the preshock (not shown), suggesting that the population is dominated by photoionization and radiative excitation from the shock.  \nThe second row of Figure 4 shows the line source function ( S l ; Equation 15). The source function varies slowly with the distance in the postshock with only a small slope, whereas in the preshock, the source function peaks at a certain distance and drops off quickly closer/further from the shock. The Planck function is plotted as dotted lines in the second row for comparison. The condition differs from LTE throughout the postshock. On the other hand, the source function is essentially the Planck  \nfunction at the outer end of the preshock, further away from the shock, as conditions approach LTE.  \nThe line-center opacity ( χ 0 ,l ) is shown in the third row of Figure 4. Similar to the population (first row), χ 0 ,l peaks in a narrow region in the postshock, whereas it is more constant in the preshock until a fast drop off. In the postshock, the opacity is higher than in the LTE cases when the flow is closer to the shock and lower than in the LTE cases when approaching the heated photosphere. In the preshock, it is higher than LTE at all distances. The resulting line-center emissivity, η 0 = χ 0 S l , is shown in the last row of Figure 4. The postshock has largely higher emissivity than the preshock. However, this does not translate to intensity.  \nFigure 5. The wavelength-dependent, two-dimensional line properties for the same model as in Fig. 4 for cos θ 0 = 0 . 5. The values in the upper row are the local properties of the lines, whereas the values in the lower row are cumulative and include the contribution from the continuum. The white dashed vertical lines mark the rest-frame line centers, and the horizontal lines mark the shock location. The red dotted lines show the shifted line-center wavelength λ ' = λ 0 ( v ( r ) cos θ 0 /c +1).  \n<!-- image -->  \nIn Figure 5, we show the wavelength-dependent line properties of the same model as in Figure 4, assuming µ ≡ 1 / cos θ 0 = 2, but the results are similar regardless of angle. The top row of Figure 5 shows the opacity χ ν,l and the emissivity η ν,l as a function of depth. The opacity and emissivity peak near the rest-frame line centers (vertical dashed lines) in the postshock, whereas a velocity shift is seen for the preshock, as expected. They are also broader and higher, suggesting that conditions in the postshock are more suitable for line formation than that in the preshock. Specifically, the higher temperatures mean that the lines are broader (higher thermal width), and the higher density means higher χ and η .  \nThe bottom row of Figure 5 shows the optical depth from the observer toward the star along the flow (left panel) and the differential intensity from each depth (right panel). The continuum opacity and emissivity are included in these calculations. Even though the opacity in the preshock is lower than in the postshock, the line optical depth is much higher in the preshock than in the postshock due to a larger physical distance. The optical depth is low throughout the postshock near the rest-frame line centers.  \nThe differential intensity dI ν = ( η l + η c ) e -µτ ν µdr (c.f., Equation 9) is shown in the lower right panel of Figure 5. Due to the high optical depth, the preshock emission near the shifted line center λ ' = λ 0 ( v ( r ) cos θ 0 /c +1) (red dotted lines) is absorbed for nearly the entire flow. Only the radiation from the two wings of each line can escape the flow and contribute to the final line emission. On the other hand, the low optical depth in the postshock allows the radiation to escape easily. Again, we note that the differential intensity in the preshock is lower than in the postshock, but its sum is higher due to the much larger thickness.  \nLastly, the Figure also shows that the emission from the preshock comes mostly from the front of the flow (i.e., closer to the observer), but most emission from the postshock comes from the back (i.e., the interface between the postshock and the heated photosphere beneath). This is where the temperature has decreased enough for C IV to be the dominant ionization stage, as shown in Figure 3. This is also where the velocity in the flow is lowest, which explains the smaller shift in velocity space from the narrow component of the line as seen in the observations (Ardila et al. 2013).  \nFigure 6. The change in radiative properties of the flow due to turbulence velocities. The units are the same as in Fig. 4.  \n<!-- image -->", '3.2. The Effects of Turbulence': 'Multidimensional simulations of accretion shocks have shown that the process can be very turbulent, depending on the ratio of gas/magnetic pressure. In the preshock, the magnetic pressure is much stronger than gas pressure, and the mass follows the field lines, so turbulence should not be dominant. In the shock and the postshock, the kinetic energy of the bulk flow is converted into heat and/or turbulence. The hydrodynamic pressure overcomes the magnetic pressure, and turbulent flows can develop (Orlando et al. 2010).  \nNevertheless, the level of turbulence in the flow is unknown, so we model it as inherent turbulence velocity. We investigate this effect on the radiative properties of the flows by varying the velocity in the Cloudy input. Figure 6 shows the line-center opacity and emissivity for the postshock and preshock. Changing the turbulent velocity does not affect the radiative properties of the postshock. In the preshock, the turbulence affects the emissivity in the entire flow and the opacity in the region where the line approaches LTE (c.f., Figure 4). This is because turbulence affects the pressure and temperature, affecting the level population as it approaches LTE further away from the radiation front.', '3.3. The Effects of Abundance': 'Observations of chemical species in protoplanetary disks suggest that the elemental abundances in the disks may be modified from the initial values due to processes such as dust trapping, grain growth, or planet formation. Variations in carbon abundance in the (inner) disk may have a non-linear effect on the profiles of the lines forming in the accretion flow/shock. Therefore, we include the abundance as another parameter in our model. Similar to the previous subsection, we show the line-center optical depth and the line emissivity for an  \nexample model of different abundance (compared to the Sun) in Figure 7. The solid lines are for models with the abundance as a parameter in the Cloudy run, whereas the dashed lines show the results when the level populations are scaled from the [C/H]=0 case by the abundance. There is no noticeable difference in the postshock for the two methods, suggesting that the formation of the C IV line in the postshock is optically thin, i.e., the curve of growth is in the linear regime. On the other hand, the preshock shows an order of magnitude or more difference in the emissivity when compared between the two methods, suggesting that the line forms in the high optical depth regime.', '3.4.1. Relative Contribution': "Figure 8 shows the contribution to the line profiles from each component of the flow structure. As seen in the left panel, the low-velocity components originate from the postshock, while the high-velocity, broad components are largely from the preshock. Nevertheless, a significant opacity in the preshock flow is still present at the wavelength of the narrow component, leading to the reduction in flux coming from the postshock. For a given spot ( θ , ϕ ), the postshock emission follows essentially a Voigt profile. In contrast, the preshock emission is strongly self-absorbed, and only the wings of the line contribute to the emission (see also Figure 5). The smooth, broad, high-velocity component is the superposition of multiple self-absorbed lines at different observed radial velocities from different ϕ . Therefore, in our model, the final line profile is sensitive to the choice of geometry.  \nThe line center of the narrow component is slightly redshifted, comparable to observations (Ardila et al. 2013). However, there is very little difference in line flux  \nFigure 7. The change in radiative properties of the flow due to different carbon abundance. The solid lines are the properties with the abundance as an input in Cloudy. The dashed lines are the values computed by scaling the level population from the results of the [C/H]=0 calculation, demonstrating the non-linear relationship between emissivity and abundance. The units are the same as in Fig. 4.  \n<!-- image -->  \nand shape from different ϕ . This is because the terms cos θ 0 cancel out each other in Equation 7 and Equation 9. The low opacity/optical depth in front of the emitting region (near the heated photosphere; Figure 5) makes the exponential term in Equation 9 insensitive to the viewing angle θ 0 .  \nThe center panel of Figure 8 shows the relative contribution to the overall line profile from the postshock and the preshock for F = 10 11 erg s -1 cm -2 . The fluxes of the narrow components are attenuated by the preshock, especially at the red side of the line. This attenuation comes from the blue side of preshock's line and the continuum. Depending on the inclination, the attenuation from the preshock may be significant, up to ∼ 50% in some cases.  \nThe total line profile of the doublet is shown in the right panel of Figure 8. As mentioned earlier, the accretion flow is split into five flows (Figure 1), with the low and moderate energy-flux flows split into the top and bottom components. To illustrate the relative importance of each flow, the filling factor for each energy flux is set to be f = 0 . 01. It is evident that the high energy flux F = 10 12 erg s -1 cm -2 contributes the most to the line profile for the same f . On the other hand, the lowest energy flow contributes the least for the same f . For the model of a real star, however, the relative contributions of the components will be different since the relative values of the filling factors are different.  \nAs noted earlier, we assume that the outer, lowerdensity parts surrounding the inner, denser part of the accretion funnel do not add much opacity and will not affect the emerging C IV profiles of the denser part. Since the shock locations depend on the flow's ram pressure and the ambient photosphere, different flows are  \nexpected to form the shock at different locations. Nevertheless, how these differences affect the C IV profiles is beyond the scope of this paper, as this requires full 3D modeling.", '3.4.2. Geometry Effects': "As mentioned earlier, the overall line shape strongly depends on the geometry. The information is encoded in the values of θ 0 , which depend on the relative values of inclination i and the co-latitude of the flows θ (c.f., Figure 1). Again, we consider the structure of the flow and the filling factors as in Section 3.4.1. For R i = 3 . 0 R ∗ , we get a co-latitude of ∼ 35 · .  \nFigure 9 shows the total line profiles of the same model for different viewing inclinations. Both the narrow and broad components shift to the blue as the inclination increases. For i < θ , the fluxes of the narrow components are high and barely change. On the other hand, the broad components show a significant reduction in flux between i = 0 · and 20 · . This is because, at low i , a wider range of ϕ is seen, and they have a similar cos θ 0 . This leads to a generally higher flux at low i as more surface area is seen. A similar cos θ 0 leads to the postshock line being in a narrow range of velocities. Therefore, the blue wings of the line in the preshock do not have much overlap over the postshock, so the emission from the postshock (which is insensitive to cos θ 0 ; Section 3.4.1) is largely unaffected. Due to high optical depth, the preshock is affected by a slight change in cos θ 0 , resulting in a significant change in fluxes.  \nFor i > θ , the flux in the preshock component continues to decrease with increasing i . The velocity shift of the preshock emission is also more prominent, such that at i = 90 · , the blue wing merges with the emission  \nFigure 8. The contribution to the line profiles of an example model with f = 0 . 01 for all three values of F , i = 60 · , v t,pre/post = 50 kms -1 , [C/H]= -1, and distance of 360 pc. The dotted lines mark the rest-frame line centers. Left: The contribution to the flux from each ϕ . The solid lines are from the preshock, and the dashed lines are from the postshock. Note that the observer is at ϕ = π/ 2. Center: The contribution of the postshock and preshock for F = 10 11 erg s -1 cm -2 . The dashed lines show the individual emission and their sum, ignoring the attenuation effect of the preshock on the postshock. The spectra plotted with the dotted lines take into account the attenuation of the postshock and the opacity of the preshock. Right: The contribution to the overall line profiles from each energy flux.  \n<!-- image -->  \ncoming from the postshock and boosts the narrow component's flux. The postshock component is unaffected by the optical depth effects, and the flux changes are due mostly to different coverage in ϕ .  \nThe right side of Figure 9 shows the changes in the line profiles for different truncation radii while other parameters are fixed. Only a slight change in the overall line fluxes is present, but the line shape remains very similar. This suggests that the shape of the profiles is more affected by the apparent shape of the emission region (controlled by the inclination) rather than the region's size or location on the star (controlled by the truncation radius).", '4. APPLYING THE MODEL TO ULLYSES ORION TARGETS': "We test whether the C IV line originates in the accretion shock by applying our model to a set of contemporaneous UV to infrared spectroscopic observations. These data were acquired through the Hubble Ultraviolet Legacy Library of Young Stars as Essential Standards Director's Discretionary program (ULLYSES; RomanDuval et al. 2020, along with community-led programs Outflows and Disks around Young Stars: Synergies for the Exploration of Ullyses Spectra (ODYSSEUS; Espaillat et al. 2022) and PENELLOPE (Manara et al. 2021). These programs observed a large sample of T Tauri stars across different mass ranges, ages, and mass accretion rates. These observations include mediumresolution FUV spectra ( HST /COS), low-resolution NUV/Optical spectra ( HST /STIS), as well as con-  \ntemporaneous ground-based high/moderate-resolution spectra from the ESO's VLT. The data are public and are accessible through the Mikulski Archive for Space Telescopes (MAST) 1 .  \nThe first set of ULLYSES stars that were studied were those in the Orion OB1 and σ Ori regions (Manara et al. 2021; Pittman et al. 2022). We selected a subset of these stars, namely CVSO 104, CVSO 107, CVSO 146, CVSO 176, CVSO 58, CVSO 90, and CVSO 109A, to be included in our study since their accretion properties have been characterized using accretion shock models (Pittman et al. 2022). The FUV spectra analyzed in this paper are obtained directly from the ULLYSES Public Data Release 6 (DR6). The COS spectra are binned at 3 pixels but are not corrected for reddening as we also fit for the A V (see Section 4.3.2).  \nThe study provides most of the required stellar and accretion parameters for our modeling, leaving only the truncation radii unknown. Optical spectra, specifically of the H α and H β lines, are needed to determine a truncation radius R i using the magnetospheric accretion model. We use contemporaneous spectra observed by the PENELLOPE program (Manara et al. 2021) and apply the magnetospheric flow model of Muzerolle et al. (2001) to determine the truncation radius of each star.", '4.1. Determining Disk Truncation Radius': "1e 13  \n1e 14  \nFigure 9. Left : The total line profile for an example model as seen at different inclinations (left) while other parameters are the same as in the right panel of Fig. 8. Right : An example of a model profile calculated at different magnetospheric truncation radii and otherwise the same accretion/geometric parameters. The dotted lines mark the rest-frame line centers.  \n<!-- image -->  \nThe truncation radius is an important input for the shock model and our line model; it is typically assumed to be 5 R ∗ (e.g., Gullbring et al. 1998; Calvet & Gullbring 1998). However, recent studies have shown evidence of this radius being less than 5 R ∗ (e.g., Espaillat et al. 2022; Thanathibodee et al. 2023). Since our model relies upon the different infall velocities among different accretion flows, better estimates of the truncation radii are needed. We note, however, that the freefall velocity only varies with the truncation radius as v s ∝ √ 1 -R ∗ /R i , so we can afford a large uncertainty in R i .  \nHere, we use the magnetospheric accretion flow model of Hartmann et al. (1994, 1998) and Muzerolle et al. (2001) to fit the hydrogen Balmer H α , H β , and H γ lines of each star. The model's parameters include the mass accretion rate ( ˙ M), truncation radius ( R i ), accretion flow width at the disk ( W r ), maximum flow temperature ( T max ), and inclination ( i ). For each star, we create a grid of models varying the mass accretion rates by a factor of a few around those found by the shock model to account for the unknown calibration of accretion rate measurements between the two models. Specifically, we vary ˙ M between 1 -5 × 10 -8 M ⊙ yr -1 for CVSO 58, CVSO 90, CVSO 104, and CVSO 107, between 1 -8 × 10 -8 M ⊙ yr -1 for CVSO 176, and between 0 . 2 -1 . 1 × 10 -8 M ⊙ yr -1 for CVSO 146. We adopt wider ranges for other parameters. Depending on the star, the grid size ranges from 6000 to 18000 models. Following the procedure outlined in Espaillat et al. (2022), we fit the resulting model profiles with the observation, taking into account the instrumental resolutions, by calculating the χ 2 . The final model parameters of a given star are  \ncalculated from the weighted average parameters from each epoch/instrument and are shown in Table 2. Appendix A shows the details of the fitting and the results for individual lines. For CVSO 109A, we adopted the flow best-fit parameters from Espaillat et al. (2022).", "4.2. Redetermining the Accretion Flow's Filling Factors": 'The filling factors of the targets determined by Pittman et al. (2022) assumed a disk truncation radius of 5 R ∗ for all but CVSO 109A. For self-consistency, we followed the same procedure outlined in that work to redetermine the filling factors of the accretion flows using the new truncation radii. In short, we calculated the Calvet & Gullbring (1998) accretion shock models for each star given its stellar parameters and R i determined according to Section 4.1. We then use an MCMC process to fit the models to the HST /STIS NUV-optical continua by varying the filling factors f h,m,l and V-band extinction A V . The re-determined mass accretion rates, filling factors of the flows, and A V are listed in Table 2. The results are slightly different from that of Pittman et al. (2022). Most notably, the accretion rates are higher in this work than in the previous work. This is to be expected since the luminosity of the accretion column L scales as L ∝ (1 -R ∗ /R i ) L acc ∝ (1 -R ∗ /R i ) ˙ M (Calvet & Gullbring 1998). Therefore, a smaller R i than the previously assumed 5 R ∗ would lead to larger ˙ M for the same observed L .', '4.3. Fitting the C IV line': 'Table 2. Results of the Accretion Flow and Shock Models', '4.3.1. Grids of Models': 'For a given star, we create a grid of models varying the pre/postshock turbulent velocities from 5 to 100 kms -1 , the carbon abundance [C/H] from -3.0 to 0.0, and the inclination ranging between those measured with the v sin i (Pittman et al. 2022) and the flow modeling (Section 4.1) with a step of 5 degrees. The geometry, which is controlled by the truncation radius (R i ) and the filling factors ( f l,m,h ), is fixed. We assume that the abundance is constant throughout the structure and that the preshock and postshock turbulence are independent. The final model profiles are convolved to the line spread function (LSF) of COS and re-binned to the same wavelength scale as the observations. Depending on the range of possible inclinations in each star, there are 10,400 to 41,600 models per grid.', '4.3.2. Fitting procedure': "The COS observations need to be corrected for the intrinsic radial velocity. Since there are no photospheric lines in the FUV spectra, we use H 2 lines to determine the radial velocity since they form in the disk (Herczeg et al. 2004; France et al. 2012; Hoadley et al. 2015; Gangi et al. 2023), and the line center should be at zero velocity compared to the star. For each line, we fit a Gaussian and determine the velocity shift of the line compared to the line center at rest. We use eight strong lines for these fits, at 1446.120 ˚ A, 1467.08 ˚ A, 1489.570 ˚ A, 1500.450 ˚ A, 1504.760 ˚ A, 1524.650 ˚ A, 1556.840 ˚ A, and 1562.410 ˚ A, and the final radial velocity is the mean of the fits. We then correct the spectra with these RVs. The fitted radial velocities for the stars are between ∼ 10 -20 kms -1 , which is comparable to the wavelength uncertainty of the COS spectrograph at ∼ 10 kms -1 .  \nThe blue side of the C IV doublet overlaps with a fluorescent H 2 emission line at 1547.34 ˚ A, which we account for in our fitting. We assume that the feature can be rep-  \nesented by a Gaussian line profile (France et al. 2012); however, line blending makes it challenging to directly measure its best-fit parameters. Instead, we first determined the best-fit Gaussian amplitude, FWHM, and central velocity of the well-separated H 2 emission line at 1489.57 ˚ A, using an interactive model fitting routine that convolves the 'intrinsic' Gaussian with the wavelength-dependent line-spread-function of the HST -COS instrument (Arulanantham et al. 2018). This feature originates from the same upper level and corresponding disk gas population as the 1547.34 ˚ A transition ([ v ' , J ' ] = [1 , 4]). Since UV-fluorescent H 2 emission lines are primarily Keplerian-broadened, with flux ratios that scale with the ratio of Einstein B coefficients, we assume that the two features share a Gaussian FWHM and central velocity. We then multiply the best-fit amplitude of the 1489.57 ˚ A line by the ratio of Einstein coefficients of the two transitions to estimate the Gaussian amplitude of the 1547.34 ˚ A emission line.  \nSince the continuum calculated from the model only includes the contribution from the shock emission (the star's photosphere is negligible), it does not represent the real observation, which may include other sources, such as the unknown chromosphere and the underlying molecular dissociation continuum (France et al. 2017, 2023; Espaillat et al. 2022). Therefore, we only fit the residual profile by subtracting the continuum in the observation and the model. We also consider the possibility that the chromosphere contributes to the C IV line emission. Ardila et al. (2013) presents the line of a number of non-accretion young stars, showing that the profiles are narrow and weak. We calculate the flux per stellar surface area around (2 ± 1) × 10 5 erg s -1 cm -2 by averaging across four stars (RECX 1, TWA 7, LKCa 17, and HBC427) using stellar parameters from Ingleby et al. (2013). Multiplying by the surface area of the stars we model in this study, we expect that the chro-  \nsphere contributes only 0 . 7 -7% to the total C IV flux we observed, thus considering it negligible.  \nGiven that the extinction at 1550 ˚ A is high, a slight change in A V may significantly change the flux level. Therefore, we also allow A V to vary between the range of those calculated from the shock modeling (Section 4.2, Table 2) and the values from fitting the optical spectra Manara et al. (2021). We expand the A V range by 0.1 magnitude on each side to account for any possible uncertainty. For model i , the fitting flux can be written as  \nF λ,i,fit = 10 -A V A ( λ ) /A ( V ) 2 . 5 F λ,i, 0 + F λ,H 2 , (18)  \nwhere F λ,H 2 is the Gaussian H 2 profile determined above. Here, we use A ( λ ) /A ( V ) from Cardelli et al. (1989) with R V = 3 . 1, which is a fine approximation for the Orion OB1b sub-association the targets belong to (Brice˜no et al. 2019). The best fit A V for each model is determined using Scipy's curve fit function. We then calculate the χ 2 for each model by  \nχ 2 i = ∑ ( F j,i,bestfit -F j,obs ) 2 σ 2 j,obs , (19)  \nwhere j indexes over the pixels within 750 km s -1 of the mean line center at 1549.488 ˚ A. The best-fit parameters for each star are the averages from the top 100 best fits, representing the top ∼ 0 . 24% to ∼ 1% of the model grid, depending on the star.", '4.4. Fitting Results': "Figure 10 shows the top 100 best fits for the sample. The model profiles roughly agree with the observations for CVSO 176, CVSO 109A, CVSO 107, and CVSO 104. In particular, the model qualitatively reproduces the narrow and broad components of the doublet's profiles and the absolute flux levels. We note that the model also reproduces the blue sides of the narrow components and the sharp drop on the red sides of the broad components.  \nSelecting 100 best fits for CVSO 176, CVSO 109A, CVSO 107, and CVSO 104, we show in Figure 11 the combined distributions of the fitted parameters The vertical lines show the median of the distribution: [C/H]=0.8, v t,pre = 35 kms -1 , and v t,post = 65 kms -1 . These indicate that the accretion flows of these four stars are carbon-poor and highly turbulent. For the preshock, the turbulent velocity is ∼ 20% or less of the free-fall velocities, with the peak around 50 km s -1 . On the other hand, the postshock seems to have a bimodal distribution in the turbulent velocity. The first peak at ∼ 30 kms -1 , is comparable to the velocity near the base of the postshock (about 90% from the shock front). The second peak near 100 km s -1 could be the artifact of the fitting where the model tries to compensate for the large width of the narrow component that it cannot reproduce by increasing the turbulence to broaden the line. This is discussed further in Section 5.2.  \nIn Figure 12, we compare the best-fit inclinations and A V with their possible ranges. There is a slight preference of inclinations toward the values determined from v sin i (Pittman et al. 2022). On the other hand, there is no clear preference of A V toward either the shock modeling or the optical spectra fitting (Manara et al. 2021).", '5.1. Geometry and Filling factors': "Analyses of a large number of CTTS by Ardila et al. (2013) showed that the shapes of the C IV doublet were diverse and did not correlate with accretion rates, line luminosity, or inclination. They also did not see any dependence of accretion rates or accretion luminosity on inclination, but the sample size was too small to be definitive. Nevertheless, their data suggested that the emitting region of the C IV doublet is inconsistent with either being localized or covering the entire star. That is, it cannot be a small spot or a very large spot. This, however, does not rule out multiple accretion spots distributed across the surface.  \nOur suggested geometry of the accretion ring is consistent with this picture. Specifically, the emission ring should be mostly seen by an observer regardless of inclination. In reality, the ring is likely not axisymmetric, but as long as the C IV emission comes from a sufficiently large range of azimuthal angles, C IV is always seen in accreting stars.  \nOne major difference between our results and Ardila et al. (2013)'s conclusion is that we found that the broad component of the C IV doublet is optically thick, similar to what was suggested by Lamzin (1998, 2003a), whereas the narrow component is optically thin. We note that, as pointed out by Lamzin (1998), the flux ratio of 2:1 between the doublet component does not necessarily mean that the emission is optically thin, but it rather reflects the electron's collisional excitation coefficient. Moreover, the self-absorption of the emission coming from an area in our model suggests that most of the local emission in the broad component is not observed, so the total line fluxes should not necessarily correlate with other properties. This is also found to be the case by Ardila et al. (2013).  \nIn our model, the broad component's shape depends greatly on the flow's geometric parameters, such as the inclination and truncation radius, as seen in Figure 9. This is because the broad component is composed of many different wings of high-optical depth lines, the shape of which also depends on the geometry. Our adapted geometry for the accretion hotspots, i.e., axisymmetric rings, may not be able to fully explain the observations of every star in our sample. Rather, the line profiles may originate from a more complex distribution of small(er) accretion hotspots on the stellar surface as described by Lamzin (2003a). Hence, we plan to explore other geometries, such as banana-shaped hotspots (Kurosawa & Romanova 2013), in future works.  \nFigure 10. The top 100 best fits. The blue lines are the observations. The red dashed lines are the Gaussian representation of the H 2 feature. The black lines are the reddened C IV models plus the Gaussian. The continua have been subtracted for the models and the observations. The models fit the observations well for stars in the top rows (CVSO 176, CVSO 109A, CVSO 107, and CVSO 104).  \n<!-- image -->  \nFigure 11. The distributions of fitting parameters from the top 100 best fits of CVSO 176, CVSO 109A, CVSO 107, and CVSO 104. The vertical lines mark the median of the distribution.  \n<!-- image -->  \nThere are also indications that, in reality, the flow of the accreted material is not just oriented towards the star as prescribed in a 1D infall model. Instead, the post-shock material hits upon the much higherdensity photosphere, and the thermodynamic pressure overcomes the magnetic confinement from the dipolar field. Ionized material then escapes sideways out of the accretion funnel. Orlando et al. (2010) see this scenario in magneto-hydrodynamic simulations and Brickhouse  \net al. (2010) invoke it to explain the O vii emission measure observed in X-rays from TW Hya, which is much larger than expected in the post-shock zone. If this scenario is true for the targets studied here, there can be additional localized, blue-shifted C IV emission or absorption (depending on geometry and temperature of the gas spilling out sideways and upwards), which is not represented in our 1D flux model. In particular, in CVSO 146, CSVO 90, and CVSO 58, there is consider-  \nFigure 12. The best fits inclination and A V compared to the values from the literature and accretion modeling. The error bars are the standard deviations from the top 100 best fits. The green shades mark stars that have good fits. The model inclinations are constrained to be between those from v sin i (Pittman et al. 2022) and from accretion flow modeling (Section 4.1). The ranges of A V are from Manara et al. (2021) and the shock modeling (Section 4.2), plus 0.1 magnitude for uncertainty.  \n<!-- image -->  \nable blue-shifted emission that is not seen in our model (Figure 10) but might be explained this way.", '5.2. Turbulence': "While the implied high turbulence velocity may be surprising given that hotspots are observed to have ∼ kG magnetic field strength (e.g., Johns-Krull et al. 2013), it is not clear if that is true for all spots at all times. In fact, observations and models have shown patterns of hotspots on the stellar surface that are more complex than what would be expected from a dipolar infall (Gregory et al. 2008; Donati et al. 2011a,b, 2024; Orlando et al. 2010; Sacco et al. 2008). Three-dimensional simulations, while not complete, have discovered many different modes of accretion with complex accretion structures (e.g., Romanova et al. 2012, 2021; Zhu et al. 2024). Therefore, it is possible that what we are seeing as turbulent flows may be the consequence of such complex systems.  \nThere is no theoretical prediction for the exact turbulence value to expect. Observationally, our model fitting strongly prefers 50 km s -1 , and this figure is not inconceivable for the preshock with a flow velocity of a few hundred kms -1 . Furthermore, slight variations of the radius where the material is loaded onto the field lines already lead to velocity variations on the 10 km s -1 scale. While the postshock's flow speed and sound speeds are  \nlower than the supposed turbulent velocity, the plasma in the region is dominated by hydrodynamic pressure, and the flow can be deflected to the side and is observed as turbulent flow as shown in 2D/3D simulations (Orlando et al. 2010; Bonito et al. 2014). In this case, the turbulent velocity can be larger than the radial flow speed or the sound speed.  \nAlternatively, it is also possible that the fitted turbulent velocity is artificially high due to the broad observed velocity profile of the narrow component compared to the calculated model at a given position on the star. In our model, the narrow component forms in the transition region between the postshock and the heated photosphere, which has very low to zero velocity. If, instead, the line forms in the upper part of the postshock, with velocity up to v freefall / 4, then the turbulent velocity does not need to be high for the model to fit the observation since the narrow component of the model profile can be broadened by the wide distribution of the observed angles. As the line forms at temperatures of × 10 5 K (c.f., Figure 3), this could indicate the postshock cools faster than (or is not as hot as) our shock model. This can be made possible if the shock lands at an angle (instead of perpendicular) to the surface, and the energy can dissipate to the side.  \nAnother possibility is that the fitted turbulence is, in fact, reflecting a range in the observed radial velocity  \ncoming, for instance, from horizontal velocity components of the flows observed at an angle. Therefore, this extra velocity component will help shift the emerging line in velocity space for accretion flows on the side of the stellar disk, making the observed line broader. This extra velocity component is not possible in our perpendicular flow model, but it is possible in slanted accretion flows.  \nBurdonov et al. (2020) conducted laboratory experiments to study how the accretion flow/accretion shock would behave in the case of a slanted flow compared to a perpendicular flow. They found that slanted flows result in extra material being ejected to the side, creating large perturbations and decreased heating. This scenario is consistent with our results. Specifically, the high turbulent velocity in our model results could indicate strong perturbation in the flow or from spreading material. Decreased heating would also cause the C IV line to form higher in the flow where the velocity is still high, making the line spread more prominent.", '5.3. Carbon Abundance': "Elemental photospheric abundances of young stars derived from their optical/infrared spectra are consistent with solar or close-to-solar values (Biazzo et al. 2011; D'Orazi et al. 2011; Spina et al. 2017). Photospheric abundances are available for 5/7 stars in our sample from APOGEE DR16. Jonsson et al. (2020) found that CVSO 104, CVSO 107, and CVSO 109 photospheric carbon abundances are consistent with slightly sub-solar values of [C/H] ∼ -0 . 2 to -0 . 4. On the other hand, CVSO 90 and CVSO 58 are more metal-poor, with [C/H] of -1 . 6 and -0 . 6, respectively. Since the model line flux relies on the postshock density, temperature structure, and shock geometry, our model simplification and assumptions may prevent us from accurately measuring carbon abundance in accretion flows. Nevertheless, the low carbon abundance required in our models to explain the C IV profiles of CVSO 104, 107, 109, and 176, as shown in Figure 11, seems to indicate that the flows' carbon abundance is even lower and, thus, substellar. These results suggest carbon depletion in the inner disk of these stars, compared to that of the whole disk or the initial cloud.  \nStudies of a neutral carbon line in the infrared (McClure 2019; McClure et al. 2020) suggest low carbon abundance in the inner disk of some young stars. These authors have found that in the ∼ 2 Myr Taurus starforming region, carbon in the inner disk can be depleted by up to a factor of 42 and by a factor of 100 in the 10 Myr TW Hya system. Our resulting carbon abundances (Figure 11) are consistent with these values, suggesting that the processes that remove or trap carbon from reaching the inner disks have already occurred in the 5 Myr Orion systems studied here.", '5.4. Resonance Scattering of the C IV Doublet': 'Due to its resonant nature, photons near the C IV doublet can interact with C IV atoms. During this interaction, an electron in the ground state is excited to the P state, leading to the absorption features of C IV, referred to as self-absorption in this study. Furthermore, when the excited electron returns to the ground state, it re-emits C IV emission photons in a different direction, a process known as resonance scattering. Therefore, considering resonance scattering instead of self-absorption enhances the strength of the C IV doublet in the other line of sight.  \nTheoretically, the flux ratio of C IV λ 1548 and λ 1550 is 2:1 due to the statistical weights of the P 3 / 2 and P 1 / 2 states. However, previous studies have reported various flux ratios of C IV doublet, ranging from 0.6 to 2.0, in objects related to white dwarfs (e.g., planetary nebulae Feibelman 1983; Vassiliadis et al. 1996, symbiotic stars Michalitsianos et al. 1988, 1992) since the flux ratio can vary in an optically thick medium for the C IV doublet. Specifically, Michalitsianos et al. 1992 determined that a flux ratio of the C IV doublet less than unity originates from a strong outflow ( > 450 kms -1 ), which is comparable to the separation of C IV λ 1548 and λ 1550 in velocity space. This variation in the flux ratio of metal doublets due to strong outflows was explored by (Chang & Gronke 2024) through the Mg II doublet near 2800 ˚ A. Consequently, if the vicinity of the shock region has an outflow velocity faster than 400 km s -1 or is optically thick for the C IV doublet, resonance scattering induces additional variations in its flux ratio and line profile. The effect of resonance scattering in the accretion shock of T Tauri stars will be explored in future works.', '6. SUMMARY AND CONCLUSIONS': "In this paper, we develop a framework for modeling the C IV doublet based on the accretion shock model of Calvet & Gullbring (1998) and the accretion flow model of Muzerolle et al. (2001). We then applied the model to a selected sample of Orion OB1 ULLYSES targets. We summarize our findings below.  \n- 1. In our model, each of the C IV lines is composed of narrow and broad components, with emission originating from the postshock and the preshock, respectively. The line is optically thin in the postshock and optically thick in the preshock.\n- 2. The line shape and line flux of the narrow component depends only on the filling factor of the shock as viewed by the observer. On the other hand, the emission from the preshock is self-absorbed, and the shape and flux are sensitive to the geometry, including the shape of the shock region and the viewing inclination.\n- 3. For the same filling factor, the line emission is dominated by the flow with the highest energy flux.  \n- 4. Our axisymmetric model can reproduce the line shapes and fluxes of 4 out of 7 stars from the ULLYSES Orion OB1 sample (CVSO 176, CVSO 109A, CVSO 107, and CVSO 104). This suggests that our model's assumption, especially the shock's geometry, could be crucial in the emerging line flux and shape. Other geometries will be explored in future works.\n- 5. Chromospheric C IV emission contributes 0 . 7-7% to the total flux in our sample.\n- 6. Both the postshock and preshock flows are turbulent. For the postshock, the high turbulent velocity could be due to extra velocity components not accounted for by the model or the postshock line emission forms at a shallower depth.\n- 7. For the four stars with good fits, we find that the carbon abundances in the flows are sub-solar and cannot be explained by the photospheric [C/H] of the stars. The results imply that the inner region of disks surrounding these four 5 Myr stars are carbon-depleted. This is consistent with inner disk carbon abundances found with different measurement methods.  \nThe authors thank the referees for their suggestions, which improved the manuscript. We thank the PENELLOPE team, particularly Carlo Manara, Laura Venuti, and Michael Siwak, for reducing the VLT data used in this work. We also thank Fred Walter for his insightful comments and suggestions.  \nThis work is supported by HST AR-16129 and benefited from discussions with the ODYSSEUS team (https: //sites.bu.edu/odysseus/); see Espaillat et al. (2022) for an overview of the ODYSSEUS survey.  \nBased on observations obtained with the NASA/ESA Hubble Space Telescope, retrieved from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute (STScI). The specific observations analyzed can be accessed via https://doi.org/110.17909/t9-jzeh-xy14. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555.  \nThis work is based on observations collected at the European Southern Observatory under ESO programmes 106.20Z8.001, 106.20Z8.002, and 106.20Z8.009; see Manara et al. (2021) for an overview of the data set.  \nThis research made use of Astropy, 2 a communitydeveloped core Python package for Astronomy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018).", 'REFERENCES': 'Ardila, D. R., Basri, G., Walter, F. M., Valenti, J. A., & Johns-Krull, C. M. 2002, ApJ, 566, 1100, doi: 10.1086/338223  \nFacilities: ESO:VLT (UVES, X-shooter, ESPRESSO)  \nSoftware: Cloudy (Ferland et al. 2017); Astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018)  \n```\nZhu, Z., Stone, J. M., & Calvet, N. 2024, MNRAS, 528,\n```  \n2883, doi: 10.1093/mnras/stad3712', 'A. FITTING BALMER LINES WITH THE MAGNETOSPHERIC FLOW MODEL': "To estimate the truncation radius of each star, we use the magnetospheric flow model of Muzerolle et al. (2001) to fit the H α , H β , and H γ lines observed with the PENELLOPE program (Manara et al. 2021). We fit each line/observation independently, and the results are shown in Table A1 and Figure A1-A2.  \nTo calculate the final parameters for each star, we first combined the fitting results from the lines observed with the same instrument on the same epoch. This is done by combining the χ 2 of a given set of parameter i as  \nχ 2 i , total = 3 ∑ j =1 χ 2 i,j χ 2 0 ,j , (A1)  \nwhere χ 2 0 is the minimum value in a given grid and j = 1 , 2 , 3 refers to H α , H β , and H γ , respectively. We then take the average of the model parameters from the top 100 best fits using L = e -χ 2 total / 2 as the weight, and the results are shown as 'Comb' in Table A1. Finally, we calculate the weighted average of the parameters from all instruments/epochs and report the results in Table 2.  \nTable A1 . Results of Accretion Flow Fitting for Individual Lines  \nTable A1 (continued)  \nTable A1 continued  \nThanathibodee et al.  \nTable A1 (continued)  \nFigure A1. The magnetospheric accretion flow model fits to individual Balmer lines. The blue lines are the observations, the grey lines are the 100 best fits, and the black lines are the average of the model profiles. The vertical/horizontal dotted lines mark the line center and the continuum, respectively. The corresponding fit parameters are shown in Table A1  \n<!-- image -->  \nmalized flux  \nNor  \nFigure A2. Same as Fig. A1.  \n<!-- image -->"}