Photophoretic Levitation and Trapping of Dust in the Inner Regions of Protoplanetary Disks
Abstract
In protoplanetary disks, the differential gravitydriven settling of dust grains with respect to gas and with respect to grains of varying sizes determines the observability of grains, and sets the conditions for grain growth and eventually planet formation. In this work we explore the effect of photophoresis on the settling of large dust grains in the inner regions of actively accreting protoplanetary disks. Photophoretic forces on dust grains result from the collision of gas molecules with differentially heated grains. We undertake onedimensional dust settling calculations to determine the equilibrium vertical distribution of dust grains in each column of the disk. In the process we introduce a new treatment of the photophoresis force which is consistent at all optical depths with the representation of the radiative intensity field in a twostream radiative transfer approximation. The levitation of large dust grains creates a photophoretic dust trap several scale heights above the midplane in the inner regions of the disk where the dissipation of accretion energy is significant. We find that differential settling of dust grains is radically altered in these regions of the disk, with large dust grains trapped in a layer below the stellar irradiation surface, in where the dust to gas mass ratio can be enhanced by a factor of a hundred for the relevant particles. The photophoretic trapping effect has a strong dependence on particle size and porosity.
1 Introduction
Dust in circumstellar or protoplanetary disks around TTauri stars contributes to an observed infrared excess in the spectral energy distribution of these sources and provides the formation material for planetary systems. The vertical distribution of dust has observational consequences and is critical for planet formation consequences in these disks (Testi et al., 2014). There are several processes at play in the inner disk that can simultaneously affect the dust distribution at a given radius. Weidenschilling (1977, 1980) considered the competition between gravity and turbulence in setting the thickness of the dust subdisk in the solar nebula as a precursor for planet formation. The beginning of what could be termed the modern treatment of the equilibrium between vertical turbulent diffusion and gravity for the dust distribution, also known as the dust settling, was in the analytical solutions for the equilibrium dust distribution given by Dubrulle et al. (1995). The work of Dullemond & Dominik (2004) solved for the time evolution of the dust distribution, to determine if specific grain sizes can be expected to settle during the lifetime of a observed disk. This previous work has thus been concerned with characterizing the appropriate form for the vertical motion of dust, and finding the equilibrium and the time evolution.
A third process, that can affect the distribution of dust is photophoresis. The form of photophoresis (Ehrenhaft, 1918) that we are concerned with in this work is the “photophoresis”, where temperature differences established across a grain by directional differences in incident light radiation, result in a force due to the interaction with the gas. Gas molecules reflecting off the hotter side of the particle pick up more energy and momentum from the hot side of the particle than the cold side, resulting in a net force on the particle directed away from the brighter incident radiation (Hidy & Brock, 1967; Rohatscheck & Zulehner, 1985; Mackowski, 1989; Beresnev et al., 1993; Rohatscheck, 1995). The reader is reminded that this force is not radiation pressure, which results from the momentum of interacting photons.
The role of photophoresis driven by direct stellar radiation in optically thin regions of the disk has been studied in a number of contexts (Krauss & Wurm, 2005; Wurm & Krauss, 2006; Krauss et al., 2007; Takeuchi & Krauss, 2008; Cordier et al., 2016; Cuello et al., 2016). Photophoretic levitation of dust above the disk surface during FU Orionis outbursts has been examined by Wurm & Haack (2009b), Wurm & Haack (2009a), in a different form and regime. This is an optically thin case where the dust is forced high enough that the particle is exposed to direct optically thin irradiation from the star and disk while the disk is in an outburst state. The resulting focus of that work is on the radial transport of large calciumaluminumrich inclusions in a layer termed an equilibrium corridor. Importantly, Wurm & Haack (2009b) recognized that the maximum height to which particles can be photophoretically levitated above the disk to is set primarily by the decreasing gas density with height, that rapidly attenuates the vertical force.
In contrast, here we follow up the suggestion in McNally & Hubbard (2015), referred to as MH15, that dust can be levitated by photophoresis below the stellar irradiation surface of the disk at moderate and high optical depths, driven solely by the disk’s own diffuse thermal radiation. We extend the treatment of photophoresis from MH15 to one consistent with the twostream radiative transfer approximation used in McClure et al. (2013, referred to as MC13). Using disk models matched to observations from MC13, that include radiative transfer, vertical hydrostatic equilibrium, and viscous heating, we calculate the photophoretic force on the particles. Then, we solve for the evolution and equilibrium of the vertical dust distribution in a protoplanetary disk under the combined influence of gravity, gas turbulence, and photophoresis, extending the method of Dullemond & Dominik (2004) to include the latter force.
In Section 2 we present the methods used in our calculations. Section 3 is devoted to introducing the novel formulation of photophoresis in the context of a twostream approximation of radiative transfer used in this work. In Section 4 we give results for a range of calculations preformed in the context of CI Tau, in Section 5 we briefly comment on V836 Tau, and finally Sections 6 and 7 contain discussion of the results and conclusions, respectively.
2 Methods
The underlying disk models are calculated with an updated version of the D’Alessio et al. (2006) code, derived from D’Alessio et al. (1998). The models used are those for CI Tau and V836 Tau from MC13. From these models, we have extracted the grid of gas density () along with the complete description of the approximate radiation field: the mean intensity () and net vertical flux of radiation ().
With this as the background, describing the gas and radiation in the disk, we perform a dust settling calculation in the style of Dullemond & Dominik (2004) in vertical columns in the disk. The conservation equation for the dust particle number volume density is
(1) 
where is the distance above the disk midplane, is the dust diffusion coefficient, and is the dust drift velocity (Dubrulle et al., 1995). Although this equation specifies the time evolution of the dust distribution, in this work we restrict our attention to steadystate solutions. The equilibrium dust distribution is calculated following Dullemond & Dominik (2004) with the modification that the advective derivative is upwinded, depending on the direction of the net force.
Dust diffusion is a consequence of the turbulent (viscous) angular momentum transport and energy dissipation in the underlying disk model. The reader should note that the literature on determining the relations, constant, and scalings has become complicated due to the use of a wide array of different definitions and assumptions. In terms of the turbulent (Shakura & Sunyaev, 1973) the radial diffusion of angular momentum is
(2) 
where is the sound speed and is the Keplerian orbital frequency around a star with mass at cylindrical radius . Note that the MC13 models assume a single constant value of for each disk. We follow the convention often used in numerical experiments that the dust vertical diffusivity is
(3) 
which serves to define as the vertical Schmidt number.
The stopping time is captured by , the Stokes number of the dust particle, being defined as
(4) 
(in the terms of Youdin & Lithwick (2007) this assumes the eddy timescale ). In this paper the particles treated are overwhlemingly in the limit by the time they have settled to their equilibrium distribution. Only when particles are initially falling from the wellmixed initial condition used in the calculations do some have at large altitudes.
Simulations of magnetorotational instability driven turbulent disk flows suggest that, for particles in the vanishing Stokes number limit, is a factor of a few larger than unity and can differ significantly from the horizontal Schmidt number (Carballido et al., 2005; Johansen & Klahr, 2005; Turner et al., 2006; Fromang & Nelson, 2009; Nelson & Gressel, 2010; Zhu et al., 2015). Moreover, as the stopping time grows to be significant, the Stokes number grows, and the vertical Schmidt number again varies. To capture both of these we combine a limiting value with the dependence on the Stokes number given by Youdin & Lithwick (2007) as
(5) 
where is the limiting value as . The constant is uncertain, but our results will have a important dependence on it. For a discussion of the physical origin of a Schmidt number different from unity in magnetorotational instability driven accretion disk turbulence the reader is referred to Fromang & Nelson (2009). In the calculations in this paper, the Stokes numbers for particles affected by photophoresis are very small, as photophoresis arises from collisions with the gas molecules, so it operates best in the vanishing Stokes number limit.
The stopping time is determined by the drag force on the dust particle. In this work the dust particles have sizes much less than the gas molecule mean free path and drift velocities much less than the gas molecule thermal velocity, so we use the expression for the drag force in the Epstein regime for spherical particles given by Woitke & Helling (2003) as
(6) 
where is the particle radius, is the thermal velocity of gas molecules, is the dust particle velocity, is the Boltzmann constant, is the gas temperature, is the mass of one , and is the average gas molecule mass in . The stopping time for a dust particle with mass is then defined by
(7) 
The final term in Equation (1) advects the dust vertically at the drift velocity . This is defined as the terminal drift speed where the vertical gravity, photophoresis, and drag forces balance
(8) 
The terminal velocity approximation is appropriate where the stopping time is small compared to the particle’s free fall time. The vertical component of the central star’s gravity on the dust particle is
(9) 
where is the stellar mass. The appropriate form for the photophoretic force is presented in the next section.
2.1 Consistency of the Modeling Approach
The approach used in this work is to use the thermal and radiative transfer structure from the MC13 models and postprocess dust settling calculations overtop of this. The MC13 models themselves already contain a prescription for the vertical distribution of small and large dust, so the combined models in this work are only strictly selfconsistent if the dust population which is distributed differently in the separate settling calculations in this work make a small contribution to the total Rosseland mean opacity. Additionally, in this work, the calculations only present the relative distribution of dust of a single particle type in a given column of the disk, and the absolute quantity and overall composition of the dust in the disk is not specified. Discussion of such further extensions to the approach, in light of the results obtained in this work, are postponed until Section 6.1.
3 Photophoresis with TwoStream Radiative Transfer
Although the midplane of the inner parts of a protoplanetary disk are optically thick, the optical depth decreases with height. If particles are photophoretically lofted high enough, they will be moved beyond the region where the optically thick approximation of MH15 strictly applies. The optically thick approximation in MH15 matches this twostream approach at high optical depths near the midplane, but underestimates and eventually gives the wrong sign for the photophoretic force in the vertical direction at high altitudes, where the vertical gradient of gas temperature reverses, while the direction of the radiative flux remains the same. As we find that dust particles do indeed reach such altitudes, we develop here the required expressions using the results of a twostream radiative transfer calculation to derive the photophoretic force. Since the underlying disk models used in this work use the twostream approximation on vertical rays for the radiation field, the photophoretic force derived in this way is exactly consistent with the radiation field in those models. Again, the medium is considered to be dilute, so that the gas molecule mean free path is much larger than the size of the particles experiencing the photophoretic force. Gas molecule mean free paths can be estimated as (Shariff, 2009)
(10) 
which, in the regions of the disks under consideration is comfortably larger than the particles’ sizes considered. The opacity of the medium is considered to be dominated by small dust grains, of size comparable to and smaller than the wavelengths of the radiation field. The analysis of the photophoretic force here is again appropriate for opaque particles with low reflectivity at the appropriate wavelengths much larger than the wavelengths of the radiation field.
In the twostream radiative transfer approximation applied to a plane parallel atmosphere, the radiation intensity field is approximated by its first and second moments, or equivalently the mean intensity and net flux in the direction as
(11) 
with being the angle away form the positive axis. Note, this expression in the optically thick limit reduces to the expression given in MH15 by their Equation (10). The linearization and solution of the problem proceed in the same way as in MH15 their Section 2.2 and 2.3, so we do not repeat those steps here. The constant terms yield the relation
(12) 
where is the coefficient for heat conduction to the gas given by MH15 in their Equation (17). This expression can be solved analytically for .
However, the analysis of the photophoretic force is made in the limit that ,
so we will proceed using that approximation, as it is generally accurate at moderate and
high optical depths.
(13) 
where is the thermal conductivity of the particle. Applying the approximation and using Equation (7) of MH15 gives the photophoretic force as
(14) 
where is the surface accommodation coefficient which we take to be as in MH15. This expression can be compared to Equation (26) of MH15, which gives the optically thick limit. Note that in Equation (14) is strictly defined as the direction of the radiation streams in the twostream approximation. The limitations of the accuracy of the photophoretic force, due to the linear approximation of the force integral over the particle surface, are similar.
As an approximation, Equation (14) is best in the optically thick and moderate optical depth regimes where Equation (11) is a good representation of the radiative intensity field. Although the twostream approximation produces the correct net radiative flux at all optical depths, Equation (14) does not capture the possible onesided nature of direct illumination which can occur at very low optical depths. We computed numerical solutions of the heat transfer problem, for a sphere illuminated from a single direction at low optical depth, and found that in that worst case Equation (14) underestimates the true force by approximately a factor of two. However, in the application presented in this paper such a situation does not occur.
The thermal conductivity of the particles plays an important role in the strength of the photophoretic force (Loesche & Wurm 2012, MH15, Cuello et al. 2016) Here, we use the same formulation as in MH15. This is based on the experiments for packed silicate aggregates of Krause et al. (2011) and has the form where , and is the grain’s volume filling factor. The volume filling factor is defined for a given sample as where is the density of the bulk sample, and is the mass density of the constituent solid silicate monomers.
The spin of dust particles can damp the photophoresis effect if it is fast enough so that a thermal gradient fails to develop across the particle. For the dust particles studied in this work, the rotation period of dust is expected to be slow compared to the thermal conduction timescale, so photophoresis is efficient. Details of this analysis are given in the appendix.
4 Results for CI Tau
Parameter  Value 

Stellar Mass ()  
Accretion Rate ()  
Disk Mass ()  
Inner Disk Wall radius, midplane ()  0.12 
CI Tau is a ‘typical’ classical T Tauri system undergoing an average amount of accretion from the disk onto the star. Vital parameters from the best fit model of MC13 as used in this analysis, are reproduced in Table 1. Figure 1 gives the density and temperature structure of the region of the disk model from MC13 studied here. Along with the density and temperature structures, two surfaces relating to the radiation are shown. The disk photosphere, as defined by the height where the optical depth to outgoing disk thermal radiation is in the Rosseland mean, and the surface where incoming stellar irradiation is scattered or absorbed at , where the optical depth to incoming stellar irradiation is in the Planck mean (D’Alessio et al., 1998). The disk is heated both by viscous heating and stellar irradiation. In this radial range, the midplane is dominantly heated by viscous dissipation, and the temperature decreases through the photosphere, from where thermal radiation can freely escape. Above this, the disk atmosphere temperature rises again as the combined effects of the low opacity of the gas to thermal radiation prevents it from cooling efficiently, and the stellar irradiation, which penetrates down to , heats the atmosphere. The relative vertical and radial gradients of the radiation field can be crudely judged by the slope of the disk photosphere, which deviates less than one degree from parallel to the midplane in the region shown. Hence, the radiation field in this region can indeed be well approximated as dominantly vertical, perpendicular to the midplane. Two dust populations are included in this model; a population of small dust which is wellmixed with the gas, and a distribution extending to large sizes which is settled to a thin layer with a dust scale height one tenth of the gas pressure scale height. The mass of settled dust is assumed to be vertically conserved, so the dust/gas ratio is decreased in the upper layers and proportionally increased in the midplane (D’Alessio et al., 2006). Due to its inclusion of viscous heating by accretion, parametrized by an viscosity , the disk temperature structure is dominated by a vertical temperature inversion at small radii. The gas temperature at moderate disk altitudes is lower than either the hot midplane below, heated by viscous dissipation, or the warm atmosphere slightly above, heated by stellar irradiation. As heat in the form of radiation can only escape vertically, the net flux of thermal radiation must be directed vertically up everywhere, resulting in a photophoresis force directed up, away from the midplane.
4.1 A Photophoretic Dust Trap
The upwards photophoretic force grows from zero at the midplane as the vertical flux of radiation increases and is rapidly cut off with increasing altitude, not by the properties of the disk radiation field, but by the rapidly decreasing gas density in the disk atmosphere. Thus, while the details of the radiation field are critical near the midplane, the gas density dominates above the disk photosphere. At the gas density steeply attenuates the upwards photophoretic force, the vertical component of the central star’s gravity is increasing linearly with altitude. This situation is demonstrated in Figure 2, where the net acceleration of an example dust particle is shown. It is positive from the midplane up to a stable equilibrium point at approximately Height/Radius, and negative at larger altitudes where vertical gravity dominates. Thus, a photophoretic trap for dust grains which experience photophoresis is formed at Height/Radius.
The drop off in the gas density in also critical to setting the dust settling height when considering grains which do not, or only weakly, experience photophoresis. However, for the very small particles the Stokes number at each height in the column is much lower than it is for the particles large enough to experience photophoresis. Hence, the large grains levitated by photophoresis are not expected to rise above the micronsized gains which provide opacity. Large grains will be shielded from direct radiation from the central star by the small grains. This situation is in strong contrast to that envisioned by Wurm & Haack (2009b) during FU Orionis outbursts, where large grains are levitated above the opacity provided by small particles and remain directly exposed to the stellar irradiation. Levitated large grains will still be in an environment where they collide with small grains, albeit at a lower frequency than at the midplane due to the lower density. Similarly, the photophoretic trap for large grains forms below the surface where the optical depth to stellar irradiation is significant (). Since the large grains are not directly exposed to the star’s light, they are thus not expected to have a direct effect on the emission of scattered light from the disk.
Figure 3 shows what we will take in this paper as the fiducial result for the vertical settling problem given by Equation 1, with solid () particles of radius and the assumption . Equilibrium solutions are shown, both with and without the photophoretic force. We plot both the dust density normalized to the global peak density of a wellmixed (constant dust/gas mass ratio) dust, and the enhancement in the local dust/gas mass ratio over wellmixed dust. The first gives the distribution of the dust in absolute terms on a linear axis scale, and the second gives the local enhancement in the dust density over an hypothetical initial wellmixed state on a logarithmic axis scale, which is useful as a figure of merit for considering stability. The result without photophoresis shows the conventional behavior of dust settling for particles with small Stokes number; that particles settle out of the very lowest density regions, and sediment below, and moreover for these small grains, the dust remains to a good approximation wellmixed up to where it has settled out.
We find for this case the local dust/gas enhancement reaches fifty to one hundred for a wide range in the inner disk. If the background wellmixed dust/gas mass ratio of dust grains which behave in this way (i.e. excluding small grains) is then the local dust/gas mass ratio would be unity in the dust trap. At levels even approaching this figure, instabilities of some manner can be expected (Youdin & Goodman, 2005; Johansen et al., 2006).
4.2 The Vertical Schmidt Number and
The vertical Schmidt number used in this work was introduced in Equation (3). Numerical experiments attempting to determine the correct value for magnetorotational instability driven turbulence have a wide range of results (Carballido et al., 2005; Johansen & Klahr, 2005; Turner et al., 2006; Fromang & Nelson, 2009; Nelson & Gressel, 2010; Zhu et al., 2015) and, as an additional complication, the unspecified instability leading to the viscosity in the MC13 models need not be magnetorotational instability. Thus, we consider here a variation of around the fiducial case of used in the previous section. Note here, how we follow a definition in relation to the radial angular momentum transport. The same calculation as in the previous section, but with is shown in Figure 4 and that with is shown in Figure 5. It is clear that the value of the Schmidt number adopted in this model has a significant impact on the trapping of dust. Generally, a Schmidt number larger than unity decreases the tendency of turbulence to return the dust to a wellmixed state at the midplane, leaving more particles in the photophoretic trap and a lower dust density at the midplane, as reflected dramatically in Figure 5. We note that this strong dependence of the dust distribution on the vertical Schmidt number is quite different from the behavior in the case neglecting photophoresis, where the dust distribution barely varies as the vertical Schmidt number in changed in this range (see Figure 4 and Figure 5, green dotted lines).
4.3 Particle Size
The outcome of photophoretic levitation varies importantly with particle size. The approximation, that incident photons deposit their energy on the side of the particle they encounter, begins to break down for particles with size of order the wavelength of the light. Here, we vary particle size, but defer the detailed treatment of the regime, where photophoresis is gradually cut off for small particle sizes, for a future work. As the peak wavelength for thermal radiation is , the smallest particles we consider have radius , a quarter the size of the fiducial case shown in Figure 6. Thus, we expect the lower limit on the size of particles photophoretically levitated to be set by the characteristic wavelength of the incident radiation field, as opposed to the relative scaling of the diffusive and photophoretic effects with particle size.
At the other end of the grain size spectrum, although the photophoretic force increases with particle radius, so does the particle mass, making dust settling more effective. Figures 7 and 8 show the effect of increasing the grain sizes from the fiducial value to and . From Figure 8, we conclude that in the case of CI Tau, the effective cutoff in particle size for the effect of photophoretic levitation in the inner disk is a fraction of a centimeter.
The size dependence of dust photophoretic levitation introduces a strong size sorting effect. In the inner disk, in addition to the enhancement in the photophoretic trap, the midplane dust is relatively depleted in particles of these intermediates sizes, due to the effect.
4.4 Particle Porosity
It is likely that silicate grains in the inner disk have a nontrivial porosity. The variation of particle porosity changes both the density of the grains and the thermal conductivity. In this section, we examine the effect of varying the particle porosity following the aggregate model used in MH15.
The most interesting regime is for large particles with significant porosity, so in Figure 9 the result for particles with radius and volume filling factor is shown, which can be compared to the result for solid () particles in shown in Figure 7. The particle trapping effect is drastically increased for those particles. For yet larger particles, with radius and shown in Figure 10 (compare with Figure 8) the effect of increasing porosity is less dramatic, although still present, as it introduces a local peak in the dust/gas enhancement in the trap. Lofting porous particles to several scale heights in this disk model presents a challenge as the collisional velocity at height can easily be above the fragmentation velocity proposed by Birnstiel et al. (2011).
We note that the results depend in some regimes on particle porosity, but varying the density of the underlying silicate monomers in the reasonable range of does not drastically affect the results considered here.
4.5 Combination of Size and Porosity
In this section, we survey a larger number of parameter values, but restrict our attention to a single radial position at . Plotting the maximum of the dust density over than range in Figure 11, in the upper panel, shows two important features. First, the local minima correspond to the minimum dust density in the photophoretic trap, falling below the increasing dust density at the midplane, as the photophoretic lofting becomes less effective. The left hand end of the curves, in particular the and cases, shows the photophoretic trapping being limited by the diffusion of particles out of the trap  as the stopping time of the small fluffy particles decreases, the maximum dust density in the trap decreases as turbulence diffuses them out. In the corresponding plot of the maximum dust/gas enhancement (Figure 11, lower panel) a local minimum is only seen in the curve, and here again it corresponds to the location of the maximum value transitioning from the dust trap to the midplane. In general, these figures show that the effect of photophoretic trapping is shifted to geometrically larger particles as the porosity increases.
4.6 Levitation/Settling Timescale
In Figure 12 the time evolution of the column of the fiducial case is given, showing that the equilibirum dust distribution is established on a timescale on the order of years. As in the case of purely gravitational settling (Dullemond & Dominik, 2004), the vertical equilibrium is established on a timescale much shorter than the accretion timescale.
5 The Case of V836 Tau
As a comparison to CI Tau, the V836 Tau disk has a much less effective photophoretic levitation effect. V386 Tau is again a classical T Tauri system, with the important characteristics of a much lower accretion rate () and a much lower implied turbulent viscosity (, MC13). Thus the disk has a lower surface density and less effective accretion heating in the inner disk. The maximum dust denstiy and dust/gas enhancement in a column at at where the midplane temperature in the model is is shown in Figure 13, and can be compared to a similar column in CI Tau shown in Figure 11. In the maximum dust density panel, the curves rise monotonically to the right, as a result of the maximum dust density always occurring at the midplane, and not in the photophoretic trap. In the lower panel showing the maximum dust/gas enhancement, the sections of the curves to the left of the minima indicate the maximum enhancement occurring in the photophoretic trap layer. In the case of V836, the photophoretic trapping effect can only be said to be mild, and even then so for the smallest and most porous grains. This indicates the role of the specific disk properties in determining the efficacy of photophoretic levitation.
6 Discussion
We have found in onedimensional dust settling calculations that in this disk model the dust density maximum occurs far above the disk midplane when the photophoretic force on the dust grains is included, in contrast to considering gravity and dust turbulent diffusion alone. Dust is lifted away from the midplane by photophoresis and stops rising when the vertical components of gravity become larger than the photophoretic force, due to the rapidly decreasing gas density of the disk atmosphere. This effect occurs only in the inner regions of the disk, where the accretion heating in the midplane and the vertical flux of energy released is significant.
The upper height of the photophoretic layer is effectively set by the decreasing gas density, or in other words the increasing Stokes number of the lofted particles. Because of this, the large particles are not found to be lofted above the very small dust grains providing the bulk of the disk opacity. Thus, although grains are lofted, the configuration proposed by Wurm & Haack (2009b) where photophoretically levitated grains enter the fully optically thin atmosphere above an FU Orionisstate disk and are directly exposed to the stellar irradiation, is not realized. This suggests that further work on that model should consider the settling and opacity caused by micronsized grains.
6.1 Directions for Further Development
Several details concerning the photophoresis of real dust particles have not been treated in this work. Importantly, the treatment of photophoresis used in this work considers only grains of an idealized type, while the underlying theory of photophoresis should be extended to include dust grain sizes comparable to the wavelength of the disk’s thermal emission, significantly nonspherical grains, and grains with very large porosity and/or nonsilicate compositions.
As discussed in Section 2.1 the models in this work can only be considered to be fully selfconsistent if the concentration of the large grains which feel photophoresis does not significantly affect the disk’s opacity. At very high concentrations this may not be justified. Indeed, when the fragmentation of large grains in the trap layer is considered there is a particular possibility for the local opacity to be modified. In this case, it would be nescessary to develop a full model for the population of dust grains of each kind, by size, composition, and porosity; solve a trial vertical settling problem for each dust species; iterate with this dust distribution for the radiativehydrostatic balance within the disk column; and then iterate, resolving the vertical settling problem and radiativehyrostatic balance problem until converged. To our knowledge, such a calculation has not been performed in the literature, even neglecting photophoresis and dust fragmentationcoagulation driven evolution. Similarly, the turbulent diffusion approximation for the dust motion is expected to break down when the dust to gas mass ratio is locally high enough (Dubrulle et al., 1995; Krijt & Ciesla, 2016). This may be achieved in the photophoretic dust trap layer if the absolute abundance of trapped grains is sufficiently high.
Localized temperature fluctuations in the disk gas can be produced by magnetohydroynamic turbulence (Hirose & Turner, 2011; McNally et al., 2014). Thermal radiation from these local hot spots locally may produce photophoresis in turn (Loesche et al., 2015; McNally & Hubbard, 2015). As these fluctuations are localized and randomized, they are expected to be similar in nature to the effect of turbulence on the diffusion of the dust, and we consider the effect on the results here to be subdominant to the uncertainty in the Schmidt number.
The disk model used here is one calibrated to match observations, but dynamical models of protoplanetary disks contain features and restrictions beyond what is included. Of greatest potential significance are the possibilities of wind driven accretion and the details of the turbulence which is modeled as a viscous . If the observed accretion rate of the system is driven by a disk wind torque (Pudritz & Norman, 1983; Gressel et al., 2015; Bai, 2016), then the inferred value of used in this paper may overestimate the turbulent energy dissipation, and hence thermal radiation flux driving levitation. The disk model used here assumes a constant value of the turbulent , whereas the relevant properties real turbulence in the disk may not have the same simple spatial dependence. Both the variation with radius and height may be different, for example due to thermodynamic or nonideal MHD effects (Bai & Goodman, 2009; Hirose & Turner, 2011). Thus, the complete understanding of dust motion driven by photophoresis in these systems ultimately demands a more comprehensive inclusion of the physics of protoplanetary disks than is currently available. At the same time, as photophoretic levitation is sensitive to these details of the interior workings of the disk, in particular the critical question of whether the accretion at small radii is driven by winds or turbulent stresses, the presence or absence of the effect may provide a way of testing disk models.
6.2 Effects on fragmentation and grain growth
The redistribution of large grains has significant implications for fragmentation and grain growth. Grain collisions driven by turbulence limit the grain size in a very porositydependent way, as solid silicate grains are much tougher than porous aggregates. In the turbulence model used in this work, the graingrain collision velocity, if dominated by the turbulent motions, is of the order of , and as the Stokes number of particles increases with the decreasing density of the disk atmosphere, and the sound speed increases with temperature, the collisional velocities are expected to grow with height (Birnstiel et al., 2009). Hence, for some porous grains the levitation could result in an increased fragmentation rate. However, the mass ratio between colliding particles will be enhanced in both the midplane and midaltitude regions by shifting the distribution of cm sized grains. With fewer collisions between grains of similar sizes, there should be less fragmentation even at relatively high velocities (Windmark et al., 2012a). This may lead to improved grain growth relative to a settlingonly case. If even a few large grains survive, they can break through the bouncing and fragmentation barriers to become seed particles for planetesimal formation (Windmark et al., 2012b). The scenario for the local concentration of the small dust population associated with the fragmentation of large grains at the midplane proposed by Krijt & Ciesla (2016) would also be altered in the presence of a photophoretic trap or even lofting effect altering the vertical distribution and collision rates of the large grains.
As in this work, we have generalized vertical settling of a single particle species to include photophoresis, the natural extension is to generalize the coupled settlingcoagulationfragmentation problem in the same way (Dullemond & Dominik, 2005). Changes to the production and distribution of small grains have the potential to alter the opacity of the disk, and have potential effects on the SED.
6.3 Impact of trap on disk dynamics
The vertical dust trap indicated by the onedimensional calculations in this work likely has important consequences when considered in higher dimensions and in fully dynamic models. As the concentration can be on the order of one hundred times wellmixed for large particles, the relative dust/gas enhancement can approach times the global dust/gas mass ratio. This enhancement suggests local dust/gas mass ratios of at least unity in the trap, or more if dust filtration in the disk has produced a higher vertically averaged dust/gas mass ratio at these inner radii. At this level, the dust momentum is no longer negligible with respect to the gas momentum, and instabilities are likely. Although as dust is lofted, the Stokes number of the particles increases due to the decreasing density. The Stokes number reached by particles described in this work is well below , and typically in the range – for the vast majority of the dust mass in the photophoretic dust trap. Thus, by the results of Carrera et al. (2015) one would not expect the layer to break up due to streaming instability. If streaming instability can be demonstrated to develop at Stokes numbers in the order of , then the photophoretic dust trap layer might provide a way of rapidly transforming severalhundredmicrometer dust directly into much larger aggregates, bypassing barriers in the coagulation growth of aggregates in the inner disk. We intend to further explore the stability of the dust trapping layer in future work.
These models are viscous disks, which at the temperatures considered are in tension with the predictions for dead zone formation from magnetohydrodynamic models (Dzyurkevich et al., 2013). They are, however, agnostic about the underlying physical mechanism giving rise to the turbulence for dust stirring and the viscosity leading to energy dissipation. An extension of this work to similar 1+1D hydrostatic disk models, which take into account the parameter dependence of dead zones, is worthwhile (Landry et al., 2013; Flock et al., 2016). If the photophoretic lofting heat source due to turbulence were to switch off due to a change in the nature of the turbulent accretion (such as the end of a gravitational instability driven deadzone bursting period), the settling of the dust layer may drive clumping through the instability proposed by Lambrechts et al. (2016).
6.4 Implications for the formation of chondrites and planets
The range of sizes for solid silicate particles which are strongly trapped is very interesting in the context of meteoritical studies. In a review of analyses of many meteorites, Friedrich et al. (2015) show the distribution of chondrule radii in H, L, and LL chondrites typically peaks near . Thus, these chondrules (reasonably approximated as solid silicate spheres) would be strongly levitated and trapped in a disk like CI Tau.
Calciumaluminumrich inclusions, if produced early in the inner part of the solar nebula, would be processed by temperature fluctuations occurring during their turbulent transport (Taillifet et al., 2013; Charnoz et al., 2015). The particle sizes, disk epoch and radii suggest photophoretic levitation has a significant effect on these histories as it will drive them away from the midplane.
Likewise, the formation of STIPs (Systems of Tightly Packed Inner Planets) has been proposed to rely on the migration and pileup of dust at small radii (Hansen & Murray, 2012; Boley et al., 2014). Indeed, even CI Tau shows radial velocity variations consistent with a hosting a hot Jupiter class planet (JohnsKrull et al., 2016). These radii overlap with those where we predict the formation of the photophoretic dust trap. Trapping grains in this layer results in a lowering the density of the same grains at the midplane, and the stability of the trapped dust layer should be explored for the potential consequences to planet formation at small radii.
7 Conclusions
In this work we present a treatment of photophoresis using a twostream radiative transfer approximation and demonstrate that the thermal radiation of the disk itself can levitate large dust grains. By combining photophoresis, turbulence, and dust settling calculations, we obtain a clearer view of the vertical distribution of large dust grains in the inner 2 of protoplanetary disks. The application of this approach to an observationally calibrated disk structure model that includes heating by both viscous dissipation and stellar irradiation processes, for “typical” T Tauri star disks, yields the following conclusions:

The addition of the photophoretic force to the dust settling calculation creates a strong vertical dust trap above the disk midplane, at altitudes under the stellar irradiation surface.

The peak dust/gas enhancement in the trap approaches times the vertically averaged dust/gas mass ratio for the respective particle type.

The degree of enhancement depends strongly on the grain size considered and the particle porosity.
The presence of a significant enhancement in the dust/gas ratio above the midplane, as well as the redistribution of large grains, may produce interesting dynamical effects and will affect grain growth and fragmentation patterns in the innermost disk. Understanding the limits of this mechanism and its implications for planet formation in this region will be the focus of future work.
Appendix A Dust Particle Spin
In the interior of the protoplanetary disk, a spherical dust grain in thermal equilibrium with the surrounding gas has a rotation rate set by the energy associated with the rotational degree of freedom of (Krügel, 2008)
(A1) 
This rotation rate has a steep dependence on the partical radius . For a temperature at the high end of the applicable range, , the rotation period of a grain with density and radius is about , whereas for a grain it is about . Under these conditions, the angular velocity kick to the smaller grain considered due to a collision with a single molecule is of the order of and the disorder time of the brownian motion is on the order of (Krügel, 2008).
For the thermal conduction time across the grain, we can follow the estimate provided by Matthews et al. (2016) of
(A2) 
Taking the heat capacity , thermal conductivity yields conduction times of for the grain and for a grain. The closest pair of rotation period and conduction times for particles in the parameter space considered here at , is for a grain with filling factor and where the rotation period is and the conduction time is . Hence, in general, we expect the temperature gradient that drives photophoresis to be sustained against the rotation of dust particles as the heat conduction times are short compared to the other timescales in the problem.
In addition to collisions with gas molecules which bring the dust spin into thermal equilibrium with the gas, dust–dust collisions can spin up dust to much higher rates. Here, we determine an orderofmagnitude criteria for the spin rate in equation (A1) to dominate the timeaveraged spin of a dust particle. Assume that the dust collision velocities are in the order of the turbulent velocities (Birnstiel et al., 2009), then the angular momentum imparted to a dust grain from the worstcase equal size collision would be in the order of
(A3) 
For a spherical grain the surface rotational velocity kick after the collision would be
(A4) 
The angular momentum which is removed from this spinning dust particle to the thermal bath provided by the gas by a subsequent collision with a gas molecule is on the order of
(A5) 
So the spinning dust will have slowed to a thermal equilibrium after a number of the order of collisions where
(A6)  
(A7) 
So, as long as the total mass of gas colliding with the spinning dust particle per unit time is much larger than the total mass of other dust particles colliding with it, then the spin will not wander far from equilibrium in a timeaveraged sense. We define this ratio of colliding gas to dust as , and it can be expressed as
(A8) 
where and is the dusttogas mass ratio locally in the disk (canonically on the order of ). So this quantity, given that , and are all small quantities, is
(A9) 
Hence, small particles can be expected to, in a timeaveraged sense, typically spin at the gas thermal equilibrium rate given by equation (A1), because they interact with a much larger mass of gas than of dust in a given time.
Footnotes
 affiliation: current address: Astronomy Unit, Queen Mary University of London, Mile End Rd, London, E1 4NS, UK
 Many similar but distinct definitions and terminology exist, in particular similar quantities have been referred to as an (effective) Prandtl number (ex: Carballido et al., 2005).
 Because of this, the additional effects considered by Loesche et al. (2016) are subdominant here.
References
 Bai, X.N. 2016, ApJ, 821, 80
 Bai, X.N., & Goodman, J. 2009, ApJ, 701, 737
 Beresnev, S., Chernyak, V., & Fomyagin, G. 1993, Physics of Fluids, 5, 2043
 Birnstiel, T., Dullemond, C. P., & Brauer, F. 2009, A&A, 503, L5
 Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2011, A&A, 525, A11
 Boley, A. C., Morris, M. A., & Ford, E. B. 2014, ApJ, 792, L27
 Carballido, A., Stone, J. M., & Pringle, J. E. 2005, MNRAS, 358, 1055
 Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43
 Charnoz, S., Aléon, J., Chaumard, N., Baillié, K., & Taillifet, E. 2015, Icarus, 252, 440
 Cordier, D., Prada Moroni, P. G., & Tognelli, E. 2016, Icarus, 268, 281
 Cuello, N., Gonzalez, J.F., & Pignatale, F. C. 2016, MNRAS, 458, 2140
 Cuzzi, J. N., Dobrovolskis, A. R., & Champney, J. M. 1993, Icarus, 106, 102
 D’Alessio, P., Calvet, N., Hartmann, L., FrancoHernández, R., & Servín, H. 2006, ApJ, 638, 314
 D’Alessio, P., Cantö, J., Calvet, N., & Lizano, S. 1998, ApJ, 500, 411
 Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237
 Dullemond, C. P., & Dominik, C. 2004, A&A, 421, 1075
 —. 2005, A&A, 434, 971
 Dzyurkevich, N., Turner, N. J., Henning, T., & Kley, W. 2013, ApJ, 765, 114
 Ehrenhaft, F. 1918, Annalen der Physik, 361, 81
 Flock, M., Fromang, S., Turner, N. J., & Benisty, M. 2016, ArXiv eprints, arXiv:1604.04601
 Friedrich, J. M., Weisberg, M. K., Ebel, D. S., et al. 2015, Chemie der Erde / Geochemistry, 75, 419
 Fromang, S., & Nelson, R. P. 2009, A&A, 496, 597
 Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84
 Hansen, B. M. S., & Murray, N. 2012, ApJ, 751, 158
 Hidy, G. M., & Brock, J. R. 1967, J. Geophys. Res., 72, 455
 Hirose, S., & Turner, N. J. 2011, ApJ, 732, L30
 Johansen, A., Henning, T., & Klahr, H. 2006, ApJ, 643, 1219
 Johansen, A., & Klahr, H. 2005, ApJ, 634, 1353
 JohnsKrull, C. M., McLane, J. N., Prato, L., et al. 2016, ArXiv eprints, arXiv:1605.07917
 Krause, M., Blum, J., Skorov, Y. V., & Trieloff, M. 2011, Icarus, 214, 286
 Krauss, O., & Wurm, G. 2005, ApJ, 630, 1088
 Krauss, O., Wurm, G., Mousis, O., et al. 2007, A&A, 462, 977
 Krijt, S., & Ciesla, F. J. 2016, ArXiv eprints, arXiv:1603.02630
 Krügel, E. 2008, An introduction to the physics of interstellar dust (Taylor and Francis)
 Lambrechts, M., Johansen, A., Capelo, H. L., Blum, J., & Bodenschatz, E. 2016, ArXiv eprints, arXiv:1604.00791
 Landry, R., DodsonRobinson, S. E., Turner, N. J., & Abram, G. 2013, ApJ, 771, 80
 Loesche, C., & Wurm, G. 2012, A&A, 545, A36
 Loesche, C., Wurm, G., Jankowski, T., & Kuepper, M. 2016, J. Aerosol Sci., 97, 22
 Loesche, C., Wurm, G., Teiser, J., et al. 2015, in LPI Contributions, Vol. 1856, 78th Annual Meeting of the Meteoritical Society, 5137
 Mackowski, D. W. 1989, Int. J. Heat and Mass Transf., 32, 843
 Matthews, L. S., Kimery, J. B., Wurm, G., et al. 2016, MNRAS, 455, 2582
 McClure, M. K., D’Alessio, P., Calvet, N., et al. 2013, ApJ, 775, 114
 McNally, C. P., & Hubbard, A. 2015, ApJ, 814, 37
 McNally, C. P., Hubbard, A., Yang, C.C., & Mac Low, M.M. 2014, ApJ, 791, 62
 Nelson, R. P., & Gressel, O. 2010, MNRAS, 409, 639
 Pudritz, R. E., & Norman, C. A. 1983, ApJ, 274, 677
 Rohatscheck, H. 1995, J. Aerosol Sci., 26, 717
 Rohatscheck, H., & Zulehner, W. 1985, J. of Colloid and Interface Sci., 108, 457
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
 Shariff, K. 2009, Annual Review of Fluid Mechanics, 41, 283
 Taillifet, E., Baillié, K., Charnoz, S., & Aléon, J. 2013, in Lunar and Planetary Inst. Technical Report, Vol. 44, Lunar and Planetary Science Conference, 2007
 Takeuchi, T., & Krauss, O. 2008, ApJ, 677, 1309
 Testi, L., Birnstiel, T., Ricci, L., et al. 2014, Protostars and Planets VI, 339
 Turner, N. J., Willacy, K., Bryden, G., & Yorke, H. W. 2006, ApJ, 639, 1218
 Weidenschilling, S. J. 1977, MNRAS, 180, 57
 —. 1980, Icarus, 44, 172
 Windmark, F., Birnstiel, T., Güttler, C., et al. 2012a, A&A, 540, A73
 Windmark, F., Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2012b, A&A, 544, L16
 Woitke, P., & Helling, C. 2003, A&A, 399, 297
 Wurm, G., & Haack, H. 2009a, in Astronomical Society of the Pacific Conference Series, Vol. 414, Cosmic Dust  Near and Far, ed. T. Henning, E. Grün, & J. Steinacker, 509
 Wurm, G., & Haack, H. 2009b, Meteoritics and Planetary Science, 44, 689
 Wurm, G., & Krauss, O. 2006, Icarus, 180, 487
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459
 Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588
 Zhu, Z., Stone, J. M., & Bai, X.N. 2015, ApJ, 801, 81