# Planetary Accretion Shocks with a Realistic Equation of State

Zhuo Chen (陈卓) Department of Astronomy, Tsinghua University
Beijing 100084, China  Xuening Bai (白雪宁) Department of Astronomy, Tsinghua
University
Beijing 100084, China  Institute for Advanced Study, Tsinghua University
Beijing 100084, China

###### Abstract

The final stage of gas giant formation involves accreting gas from the parent
protoplanetary disk. In general, the infalling gas likely approaches a free-
fall velocity, creating an accretion shock, leading to strong shock heating
and radiation. We investigate the kinematics and energetics of such accretion
shocks using 1D radiation hydrodynamic simulations. Our simulations feature
the first self-consistent treatment of hydrogen dissociation and ionization,
radiation transport, and realistic grey opacity. By exploring a broad range of
giant planet masses (0.1-3MJ) and accretion rates
($10^{-3}$-$10^{-2}$M${}_{\oplus}\cdot\rm{yr}^{-1}$), we focus on global shock
efficiency and the final entropy of the accreted gas. We find that radiation
from the accretion shock can fully disassociate the molecular hydrogen of the
incoming gas when the shock luminosity is above a critical luminosity.
Meanwhile, the post-shock entropy generally fall into “cold” ($\lesssim
12k_{\rm{B}}/m_{\ce{H}}$) and “hot” ($\gtrsim 16k_{\rm{B}}/m_{\ce{H}}$) groups
which depends on the extent of the endothermic process of H2 dissociation.
While 2D or 3D simulations are needed for more realistic understandings of the
accretion process, this distinction likely carries over and sheds light on the
interpretation of young direct imaging planets.

Accretion (14); Hydrodynamical simulations (767); Planet formation (1241)

††software: Guangqi (Chen & Bai, in prep), Matplotlib (Hunter, 2007), Petsc
(Balay et al., 1997, 2019)

##  1 Introduction

A major paradigm shift in our understandings of protoplanetary disks (PPD)
over the past few years is the ubiquity of disk substructures found in high-
resolution observations of protoplanetary disks (e.g., van der Marel et al.,
2013; Dipierro et al., 2015; Pérez et al., 2016; Andrews et al., 2018; Long et
al., 2018; Avenhaus et al., 2018). As a common interpretation, these
substructures are considered as the outcome of planetary-mass companions
interacting with the disk, which generally requires the formation of giant
planets which can naturally open gaps (e.g., Goldreich & Tremaine, 1979; Lin &
Papaloizou, 1986; Goodman & Rafikov, 2001; Dong et al., 2017), create
asymmetries (e.g., de Val-Borro et al., 2007; Zhu & Stone, 2014), or drive
spiral patterns (e.g., Dong et al., 2015; Bae & Zhu, 2018; Bae et al., 2021).
Despite the expectation of multiple planets carving substructures in dozens of
disks, including evidence from kinematic signatures in a few systems (e.g.,
Pinte et al., 2018, 2020; Teague et al., 2019), so far the only bona fide
detection has been the two giant planets in the central cavity of the PDS 70
disk (Keppler et al., 2018; Wagner et al., 2018; Haffert et al., 2019). The
two planets are found to be accreting from the parent PPD as the pathway to
build up their masses, showing H$\alpha$ emission characteristic of planetary
accretion shocks (Aoyama et al., 2018; Aoyama & Ikoma, 2019; Thanathibodee et
al., 2019; Szulágyi & Ercolano, 2020).

The ability to detect young giant planets, and direct imaging of gas giants in
general after disk dispersal, crucially depends on its initial conditions,
particularly the initial entropy that largely sets the subsequent evolution of
its luminosity and temperature as main observables. The initial entropy is
closely related to the amount of the specific energy (energy per unit mass) of
the accreted gas. In general, the smaller the fraction of specific energy
retained in the gas, the lower the initial entropy. Correspondingly, models of
giant planet evolution are classified into “cold-”, “warm-”, and “hot-start”
models that primarily differ in their initial entropy (Marley et al., 2007;
Spiegel & Burrows, 2012). The differences among these models are most
pronounced over the first to a few tens of Myrs after giant planet formation
(Berardo et al., 2017), and hence crucially affect the inference of detection
limits and estimation of planet masses.

In the standard theory of giant planet formation by core accretion (Pollack et
al., 1996), the bulk of giant planet mass is built up by accreting gas from
PPD in a runaway manner. Upon becoming sufficiently massive, the planet opens
a gap and become “detached” from the disk, forming a circumplanetary disk
(CPD, e.g., Ayliffe & Bate, 2009; Tanigawa et al., 2012; Gressel et al.,
2013). How gas accretion proceeds through a CPD is unclear, but existing
studies suggest that most gas fall into the circumplanetary region from high
latitudes from the meridional flow (e.g., Szulágyi et al., 2014; Szulágyi &
Mordasini, 2017; Fung et al., 2019). Planetary accretion may proceed through
an accretion shock either from this infalling material, or via magnetospheric
accretion from the CPD. This accretion shock is also the source of the
H$\alpha$ emission observed from the PDS 70 system.

Detailed models of planetary accretion shocks is of great significance to
understand the initial conditions of giant planets. With the accreting
material approaching at near free-fall velocity, such models are usually in
1D, focusing on energy transport through the shock and the accreting column
(Marleau et al., 2017, 2019). In this letter, we show that a realistic
equation of state (EoS), that incorporates the dissociation of H2, plays a
highly significant role in setting the post-shock entropy of the accreting gas
that will further affect the initial conditions of giant planet evolutionary
models.

This letter is organized as follows. Section 2 describes the setup and
physical ingredients of our 1D radiation hydrodynamic simulations. Major
results are presented in Section 3, and are summarized with further discussion
in Section 4.

##  2 Physical model

In this section, we describe the setup for our accreting gas giants
simulations, highlighting the combination of radiation hydrodynamics with a
general EoS that incorporates the dissociation of hydrogen molecules.

###  2.1 Governing equations

We solve radiation hydrodynamic equations in 1D spherical geometry along the
radial ($r$) direction with a general EoS, for gas accretion towards a planet
with mass $M_{p}$. The governing equations are

$\displaystyle\partialderivative{\rho}{t}+\frac{1}{r^{2}}\partialderivative{r}(r^{2}\rho
v)$ $\displaystyle=$ $\displaystyle 0\ ,$ (1)
$\displaystyle\partialderivative{\rho
u}{t}+\frac{1}{r^{2}}\partialderivative{r}(r^{2}\rho v^{2})$ $\displaystyle=$
$\displaystyle-\partialderivative{p}{r}-\frac{\rho GM_{p}}{r^{2}}\ ,$ (2)
$\displaystyle\partialderivative{E}{t}+\frac{1}{r^{2}}\partialderivative{r}[r^{2}(E+p)v]$
$\displaystyle=$ $\displaystyle\mathbb{G}-\frac{\rho GM_{p}}{r^{2}}v\ ,$ (3)
$\displaystyle\partialderivative{E_{r}}{t}+\frac{1}{r^{2}}\partialderivative{r}(r^{2}F_{r})$
$\displaystyle=$ $\displaystyle-\mathbb{G}\ ,$ (4)

where $\rho$, $v$, $p$, $E$, and $G$ are gas density, radial velocity,
pressure, total energy density, and gravitational constant, with

$\displaystyle p$ $\displaystyle=$
$\displaystyle\sum_{i}n_{i}k_{\rm{B}}T_{\rm{g}}\ ,$ (5) $\displaystyle E$
$\displaystyle=$ $\displaystyle\rho(\epsilon_{\rm
kin}+\epsilon)=\rho(v^{2}/2+\epsilon)\ ,$ (6)

where $n_{i}$ is the number density of the $i$th species (to be specified
later), $k_{\rm{B}}$, $T_{\rm{g}}$, $\epsilon_{\rm kin}$,
$\epsilon(\rho,T_{\rm{g}})$ are Boltzmann constant, gas temperature, specific
kinetic energy and specific internal energy. Radiation energy density and
energy flux are denoted by $E_{r}$ and $F_{r}$, and $\mathbb{G}$ is the rate
of energy exchange between radiation and matter, given by

$\displaystyle\mathbb{G}$ $\displaystyle=$ $\displaystyle\kappa_{\rm{P}}\rho
c(E_{r}-a_{\rm{R}}T_{\rm{g}}^{4})\ ,$ (7)

where $\kappa_{\rm{P}}$, $a_{\rm{R}}$ and $c$ are Planck opacity, radiation
constant and speed of the light. For future convenience, radiation temperature
$T_{\rm rad}$ is defined by $E_{r}=a_{\rm{R}}T_{\rm rad}^{4}$.

The radiation sub-system is closed using the flux-limited diffusion (FLD)
approximation, which relates $E_{r}$ and $F_{r}$ by

$\displaystyle F_{r}$ $\displaystyle=$
$\displaystyle-\frac{c\lambda(R)}{\kappa_{\rm{R}}\rho}\partialderivative{E_{r}}{r}\
,$ (8) $\displaystyle R$ $\displaystyle=$ $\displaystyle\frac{|\partial
E_{r}/\partial r|}{\kappa_{\rm{R}}\rho E_{r}}\ ,$ (9)

where $\kappa_{\rm{R}}$ is the Rosseland mean opacity. In this paper, we adopt
the flux limiter $\lambda(R)$ described in Levermore & Pomraning (1981)

$\lambda(R)=\frac{2+R}{6+3R+R^{2}}\ .$ (10)

###  2.2 Equation of state and opacity

Our typical simulations encompass a temperature range between $\sim
10^{2-4}$K, where $\ce{H2}$ can be dissociated and eventually ionized, and
standard ideal gas EoS becomes questionable. As a first study, we incorporate
such physics but make a simplified assumption of chemical local-thermal
equilibrium (LTE) of H2, H, and H+. The abundance of these species can be
obtained analytically from the Saha equations at runtime according to Appendix
C of Chen et al. (2019). The EoS is given by
$\epsilon=\epsilon(\rho,T_{\rm{g}})$ in an analytical from, which avoids the
use of a tabulated EoS and improves the efficiency and accuracy in our
simulations. In particular, we consider a hydrogen mass fraction of $X=0.74$,
assuming remaining mass in helium.

One important quantity we compute in this work is the gas entropy. With the
general EoS, it is given by

$s=\sum_{i}\frac{k_{\rm{B}}n_{i}}{\rho}\left(1+\frac{d\ln Z_{i}}{d\ln
T}-\ln\frac{n_{i}}{Z_{i}}\right),$ (11)

where $Z_{i}$ is the partition function for the $i$th species, given in
Appendix A.

Opacity is a crucial physical component in our model. We adopt the same
opacity tables as in Marleau et al. (2019), combining the gas opacity table of
Malygin et al. (2014) that dominates over $1500$K, and the dust opacity table
of Semenov et al. (2003) that dominates below 1100-1200K, depending on the
density. For temperatures in between, the maximum of the value of the two
tables is taken.

Note that the opacity tables assume gas and radiation temperatures are the
same. In reality, as we will see, $T_{\rm{rad}}$ and $T_{\rm{g}}$ can be
different at the Zel’dovich spike and H2 dissociation region. Here we use the
radiation temperature to obtain opacity from lookup tables, bearing in mind
the caveat which can be improved in future works.

Table 1: Simulation parameters and results.

ID | $M_{p}$ | $\dot{M}_{p}$ | $r_{\rm{in}}$ | $p_{\rm{in}}$ | $r_{p}$ | $s_{\rm{ps}}^{l8}$ | $s_{\rm{ps}}^{l7}$ | $\eta^{\rm{phy}}$ | $\chi_{\ce{H},\rm{pre}}$ | $\chi_{\ce{H},\rm{ps}}$ | $\mathcal{M}_{\rm{pre}}$ | $T_{g,\rm{pre}}$ | $T_{g,\rm{ps}}$ | $p_{\rm{ps}}$ | $\frac{\dot{M}_{p}}{{\dot{M}_{\rm{crit}}}}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
| (MJ) | $\big{(}\frac{M_{\oplus}}{\rm{yr}}\big{)}$ | (RJ) | (bar) | (RJ) | $\big{(}\frac{k_{\rm{B}}}{m_{\ce{H}}}\big{)}$ | $\big{(}\frac{k_{\rm{B}}}{m_{\ce{H}}}\big{)}$ | % | % | % |  | (K) | (K) | ($\frac{\rm{bar}}{10^{3}}$) |
1 | 0.1 | $10^{-2}$ | 1.0 | 2 | 1.284 | 11.42 | 11.42 | 82.23 | 16.67 | 5.18 | 5.24 | 1865 | 1920 | 2.850 | 0.56
2 | 0.3 | $10^{-2}$ | 1.1 | 5 | 1.294 | 17.04 | 17.06 | 64.90 | 98.73 | 71.55 | 6.00 | 2458 | 2474 | 5.077 | 1.28
3 | 1 | $10^{-2}$ | 1.5 | 10 | 1.679 | 19.21 | 19.20 | 81.59 | 99.93 | 95.72 | 8.39 | 2685 | 2754 | 4.987 | 2.46
4 | 1 | $10^{-2}$ | 1.7 | 10 | 1.865 | 18.25 | 18.20 | 82.34 | 99.72 | 83.39 | 8.61 | 2495 | 2526 | 3.765 | 1.81
5 | 3 | $10^{-2}$ | 1.7 | 10 | 1.780 | 19.61 | 19.61 | 92.81 | 100 | 99.92 | 12.57 | 3358 | 3465 | 7.843 | 7.05
6 | 3 | $10^{-2}$ | 1.9 | 10 | 1.998 | 19.68 | 19.67 | 92.12 | 100 | 99.74 | 12.37 | 3092 | 3201 | 5.791 | 4.95
7 | 0.1 | $10^{-3}$ | 1.0 | 2 | 1.266 | 11.45 | 11.45 | 91.85 | 0 | 0 | 6.06 | 1290 | 1299 | 0.304 | 0.06
8 | 0.3 | $10^{-3}$ | 1.1 | 5 | 1.198 | 11.24 | 11.23 | 97.02 | 1.24 | 0.11 | 11.82 | 1456 | 1468 | 0.696 | 0.24
9 | 1 | $10^{-3}$ | 1.5 | 10 | 1.568 | 11.49 | 11.48 | 98.52 | 11.76 | 1.21 | 18.04 | 1612 | 1667 | 0.702 | 0.37
10 | 1 | $10^{-3}$ | 1.7 | 10 | 1.785 | 11.47 | 11.48 | 98.57 | 4.43 | 0.35 | 17.83 | 1511 | 1536 | 0.489 | 0.25
11 | 3 | $10^{-3}$ | 1.7 | 10 | 1.734 | 13.16 | 13.21 | 98.34 | 71.71 | 20.41 | 23.12 | 1831 | 2040 | 1.130 | 0.81
12 | 3 | $10^{-3}$ | 1.9 | 10 | 1.940 | 12.18 | 12.21 | 98.92 | 48.3 | 8.18 | 24.02 | 1740 | 1889 | 0.817 | 0.58

Note. — From column 1 to column 16: (1) the model ID, (2) planet mass, (3)
accretion rate, (4) inner boundary radius, (5) inner boundary pressure, (6)
radius where $p_{\rm{ram}}=\rho v^{2}=p$, (7) post-shock entropy with 8 SMR
levels, (8) post-shock entropy with 7 SMR levels, (9) global shock efficiency,
defined in Equation 18, (10) the number fraction of H at the pre-shock radius,
(11) the number fraction of H at the post-shock radius, (12) pre-shock gas
Mach number, (13) pre-shock gas temperature, (14) post-shock gas temperature,
(15) post-shock pressure, (16) and the actual accretion rate divided by the
critical accretion rate.

###  2.3 Simulation setup

We solve the numerical problem with Guangqi (Chen & Bai, in prep), a new 1D
radiation-hydrodynamic code with adaptive- and static-mesh-refinement
(AMR/SMR). It employs the FLD approximation for grey radiative transfer, which
is solved implicitly and is self-consistently coupled with a general EoS. In
this problem, the shock and the planetary atmosphere require high resolution.
They are located at the bottom of the computational domain and we employ SMR
to properly resolve the Zel’dovich spike behind the shock (Section 3.3).

Our simulation domain spans between $[r_{\rm in},r_{\rm max}]$, where $r_{\rm
max}$ is fixed to 20RJ (RJ is Jupiter radius), using an uniform grid with 256
cells at base level. We use 7 or 8 levels of mesh refinement with each level
doubling the resolution of the parent level. Therefore, the finest cell has a
length of approximately 43.6 km or 21.8 km ($r_{\rm{in}}$ varies from 1 to
1.9RJ). In comparison, the finest cell in Marleau et al. (2019) is 35.7 km.

Gas is injected from the outer boundary assuming free-fall velocity, and gas
density at the outer boundary is parameterized by the accretion rate
$\dot{M}_{p}$, given by

$\displaystyle v_{\rm{inj}}$ $\displaystyle=$
$\displaystyle-\sqrt{\frac{2GM_{p}}{r_{\rm{max}}}}\ ,$ (12)
$\displaystyle\rho_{\rm{inj}}$ $\displaystyle=$
$\displaystyle-\frac{\dot{M}_{p}}{4\pi f_{c}r_{\rm{max}}^{2}v_{\rm{inj}}}\ ,$
(13)

for spherical geometry, where $f_{c}$ is the covering fraction of the
accretion flow over planetary surface, and we take $f_{c}=1$ in this work. Gas
temperature in the outer boundary is linearly extrapolated at runtime, and
radiation energy is set according to $\partial(r^{2}E_{r})/\partial r=0$.

For simplicity, we set the initial condition to be
$\rho_{\rm{init}}=\rho_{\rm{inj}}$, $v_{\rm{init}}=v_{\rm{inj}}$, and
$T_{\rm{g,init}}=T_{\rm{rad,init}}=100$ K everywhere. The particular form of
initial condition does not affect the steady state solution that we look for
after running the simulations for many free-fall timescales.

Setting inner boundary conditions requires some care. A straightforward
reflecting boundary condition would gradually build up mass and pressure near
the inner boundary, and the system hardly achieves a steady state. In reality
(especially with $f_{c}<1$), the post-shock gas flow likely spreads out,
eventually maintaining equilibrium with atmospheric pressure. This motivates
us to choose a fixed-state inner boundary condition set by
$[\rho_{\rm{in}},p_{\rm{in}},T_{\rm{g,in}}]$ in gas variables. However,
$\rho_{\rm{in}}$, $p_{\rm{in}}$, and $T_{\rm{g,in}}$ are unknown beforehand.
In practice, we choose $p_{\rm in}$ as representative atmospheric pressure at
planet surface for each simulation (to be specified in the next subsection).
We first set the inner hydro boundary condition to be non-penetrating and the
inner radiation boundary condition to be zero gradient. As the gas falls onto
the planet, we keep monitoring the increase of gas pressure at the inner
boundary. Once it reaches $p_{\rm{in}}$, we record the $\rho$ and $T_{\rm{g}}$
of the innermost cell as $\rho_{\rm{in}}$, and $T_{\rm{g,in}}$, and impose the
aforementioned fixed-state inner boundary condition. The radiation boundary
condition is still zero gradient.

Figure 1: The shock profiles of a Jupiter mass planet accreting at high (model
4, left) and low (model 10, right). From top to bottom, the first panel:
velocity and density profiles in red and black lines. Second panel:
$\kappa_{\rm{R}}$ and $\kappa_{\rm{P}}$ profiles in the red and blue, and
$L_{r}$ profile in the black line. The value of $L_{r}$ at the post-shock
radius is marked with a black circle. Third panel: $\chi_{\ce{H2}}$ and
$\chi_{\ce{H+}}$ profiles in red and blue lines. SMR levels are shown in
black. The pink color covers the region of H2 dissociation. Fourth panel:
$T_{\rm{g}}$ and $T_{\rm{rad}}$ profiles in red and blue lines, with entropy
profile $s$ shown in black. Hydrogen dissociation is an endothermic process
and marked in orange. Fifth panel: $f_{\rm{r}}$ defined in Equation 15 and
$\lambda_{\rm{mfp}}/r$ profiles in red and black lines. The radiative zone is
marked in grey.

###  2.4 Model selection

To sample a broad range of parameter space, we choose $M_{p}=[0.1,0.3,1,3]$MJ
to study the accretion of super Neptune, Saturn mass, Jupiter mass and more
massive planets. We choose
$\dot{M}_{p}=[10^{-3},10^{-2}]$M${}_{\oplus}\cdot\rm{yr}^{-1}$ to represent
the low and high accretion rates.

Since we do not model the internal structure of the planet (e.g., Mordasini et
al., 2012), which itself depends on the accretion history and shock
properties, we cannot give a self-consistent prescription of $r_{\rm in}$ and
$p_{\rm in}$ of the forming gas giants. Therefore, we treat $r_{\rm{in}}$ as
another parameter and adopt the results that newly formed gas giants with
$M_{p}\gtrsim$MJ are expected to have radii of $\sim 1.4-4$RJ (Spiegel &
Burrows, 2012). In our simulations, the accretion shock (and the Zel’dovich
spike) forms at $r_{p}>r_{\rm{in}}$, where the ram pressure ($\rho v_{r}^{2}$)
of the infalling gas equals to the pressure of the atmosphere. Note that the
free energy available in the shock scales as $M_{p}/r_{p}$ and a larger
$r_{p}$ would largely be equivalent to reducing planet mass (Section 3.3 for
more discussion). Our choice of $p_{\rm in}$ is also somewhat arbitrary,
ranging from $2-10$ bars as $M_{p}$ increases, but we have verified that the
results are insensitive to this choice. Using a larger $p_{\rm in}$ will lead
to a larger $r_{p}$ as the planetary atmosphere becomes thicker.

We carry out a total of 12 runs, varying $M_{p}$, $\dot{M}_{p}$ and $r_{\rm
in}$. The specific run parameters, as well as their major diagnostic
properties, are listed in Table 1, and they are labeled as models 1 to 12.

##  3 Simulation results

In this section, we choose two representative simulations, model 4 and 10,
corresponding to accretion onto a Jupiter mass planet with high and low
accretion rates. Results from the steady state solutions are shown and
analyzed in detail, paying special attention to the post-shock entropy111More
rigorously speaking, we refer to “post-shock” as the immediate downstream of
the Zel’dovich spike.

###  3.1 Overview of simulation results

Figure 1 shows the steady state solution of models the two models. In both
models, the gas largely free-falls onto the planet as shown in the first
panel. The infalling gas is stopped by the planetary atmosphere at $r_{p}$,
forming the accretion shock. Gas can be substantially heated to high
temperatures at $r_{\rm p}$ by the shock, followed by a rapid fall-off in
temperature, known as the Zel’dovich spike. At the Zel’dovich spike,
$T_{\rm{g}}>T_{\rm{rad}}$ (see the zoom in plot of the fourth panel), the gas
internal energy is quickly converted to the radiation energy, which is then
transported outward, seen as the emergence of high radiative energy flux in
the second panel from top. The shock is resolved by 8 SMR levels around
$r_{p}$ as shown in the third panel. We confirm that our solution is converged
in the sense that the post-shock entropy $s_{\rm{ps}}$ (to be defined later)
does not differ by much if we use 7 SMR levels.

From the large to small radii, we can approximately divide our simulation
domain into three zones: the pre-shock upstream where gas free falls; the
Zel’dovich spike where kinetic energy of the infalling gas is converted to
heat and radiation, and a radiative zone in the downstream region considered
to be a part of the giant planet upper atmosphere. The three zones are
separated by a pre-shock radius $r_{\rm{pre}}$ and a post-shock radius
$r_{\rm{ps}}$. We define $r_{\rm{pre}}$ as the radius with the maximum
infalling speed and $r_{\rm{ps}}$ as the radius where
$T_{\rm{g}}=T_{\rm{rad}}$ in the downstream of the Zel’dovich spike.

For future convenience, we define radiation mean free path as

$\lambda_{\rm{mfp}}=(\kappa_{\rm{R}}\rho)^{-1}\ ,$ (14)

and we compare $\lambda_{\rm mfp}$ to $r$ in the bottom panel of Figure 1. The
interior of the planet is optically thick with $\lambda_{\rm mfp}/r\ll 1$,
while the pre-shock region is typically optically thin.222More rigorously, one
may use the $R$ factor (Equation 9) to quantify how close the system is to the
diffusive or free-streaming limit, and the result is qualitatively the same.
We further define

$\displaystyle f_{\rm{r}}$ $\displaystyle=$ $\displaystyle F_{r}/(cE_{r}),$
(15) $\displaystyle L_{r}$ $\displaystyle=$ $\displaystyle 4\pi r^{2}F_{r}\ ,$
(16)

where $f_{\rm{r}}$ is the reduced radiation flux that characterizes how close
radiation transport is to the free-streaming limit, and $L_{r}$ is the outward
radiation flux, i.e., luminosity.

Exterior to the planet, the infalling gas undergoes compression and hence
adiabatic heating. Additional heating results from the absorption of the
outgoing radiation originating from the shock region, maintaining $T_{\rm
rad}\sim T_{\rm g}$ (Marleau et al., 2019). When the infalling gas is heated
to $T_{\rm{g,dis}}\approx 2000$K, H2 starts to dissociate. Let us define the
number fraction of hydrogen species as,

$\begin{split}\chi_{\ce{H+}}&=\frac{n_{\ce{H+}}m_{\ce{H}}}{\rho X}\
,\quad\chi_{\ce{H2}}=\frac{2n_{\ce{H2}}m_{\ce{H}}}{\rho X}\ ,\\\
\chi_{\ce{H}}&=\frac{n_{\ce{H}}m_{\ce{H}}}{\rho
X}=1-\chi_{\ce{H+}}-\chi_{\ce{H2}}\ .\end{split}$ (17)

In the third panels, the red and blue lines show the profiles of
$\chi_{\ce{H2}}$ and $\chi_{\ce{H+}}$, while $\chi_{\ce{H}}$ can be easily
deduced from Equation (17). The temperatures in these two cases are
insufficient to ionize hydrogen. For the high accretion model, H2 is almost
fully dissociated ahead of the shock. For the low accretion model, in
contrast, the infalling gas is almost molecular.

It is worth noticing that hydrogen is not ionized at the shock in model 4 but
is partially ionized in model 10. This is because the cooling strength in the
shock zone is determined by $\kappa_{\rm{P}}(a_{\rm{R}}T_{\rm{g}}^{4}-E_{r})$
and $\kappa_{\rm{P}}$ is a strong function of both density and temperature. A
high density (from higher accretion rate) at the shock leads to a large
$\kappa_{\rm{P}}$, thus a relatively low $T_{\rm{g}}$ (at the Zel’dovich
spike) is sufficient. Conversely, a low density at the shock would result in a
small $\kappa_{\rm{P}}$, and hence higher gas temperature that may ionize the
hydrogen.

####  3.1.1 H2 dissociation in the pre-shock region

By comparing the two models, we see that radiation temperature closely follows
gas temperature in the pre-shock region in general, except when H2 gets
dissociated. In model 4, the endothermic process increase the gas’s heat
capacity, leading to $T_{\rm{rad}}>T_{\rm{g}}$ in this region and hence more
radiation is absorbed by the infalling gas. This is accompanied by a drop in
$L_{r}$ between 2.5-4.7 RJ. In model 10, on the other hand, $L_{r}$ is largely
flat in the pre-shock region.

The dissociation of H2 and absorption of accretion luminosity in the pre-shock
has important consequences. The black line in the fourth panel shows the
entropy profile $s$. When $T_{\rm{g}}<T_{\rm{g,dis}}$, the gas radiates away
energy as it is compressed, resulting in a slowly decreasing entropy profile
as the gas falls inward. In model 4, upon $T_{\rm{g}}$ reaching
$T_{\rm{g,dis}}$, the dissociation of H2 and subsequent energy absorption
radiation drives the entropy to rapidly increase from 4.7RJ to 2.5RJ by almost
10$k_{\rm{B}}/m_{\ce{H}}$. This substantial increase in $s$ is not present in
model 10 with only a small fraction of H2 dissociated.

We note that at the shock, the gas internal energy is quickly converted to the
radiation energy, leading to a rapid decrease in entropy. The post-shock
entropy keeps decreasing towards the planet interior accompanied by an
outgoing radiation flux. Since there is no radiation flux at the inner
boundary, here we name the radiation flux at $r_{\rm{ps}}$ the internal
luminosity. The internal luminosity of the planets are marked with black
circles in the second panels.

We have also examined simulation results of other models. Generally, results
from models 2,3,5,6,11 and 12 are similar to model 4 because a large fraction
of H2 dissociates before the infalling gas hits the planetary atmosphere. They
differ in the radii where H2 dissociation takes place. Results from models
1,7,8 and 9, on the other hand, are similar to model 10, as they all show zero
to very low degree of H2 dissociation in the pre-shock region. They are also
very similar to the solutions of perfect gas EoS with $\gamma>4/3$ (Marleau et
al., 2017, 2019) as there is no dissociation/ionization except at the shock.

###  3.2 Global shock efficiency

Global shock efficiency quantifies the fractional energy that leaves the
system, defined as (Marleau et al., 2019)

$\eta^{\rm{phy}}=\frac{\dot{E}(r_{\rm{max}})-\dot{E}(r_{\rm{ps}})}{\dot{E}(r_{\rm{max}})},$
(18)

where $\dot{E}(r)=-\dot{M}_{p}(\epsilon+\epsilon_{\rm{kin}}+p/\rho+\phi)$ is
the total energy flux measured at $r$. A derivation of $\eta^{\rm{phy}}$ can
be found in Appendix B.

In long-term planetary evolution, shock efficiency largely determines the
initial entropy of the planet and affects the mass-radius relation of the
planet. Often, it was treated as a free parameter (Spiegel & Burrows, 2012),
and hence the outcome of the calculations somewhat hinges on this parameter.
Recently, Marleau et al. (2017, 2019) studied the global shock efficiencies
with perfect gas EoS, and conclude that $\eta^{\rm{phy}}\geqslant 97\%$ in
their models with $1.3M_{J}\leqslant M_{p}\leqslant 10M_{J}$.

In Figure 2, we can see that $64.9\%\leqslant\eta^{\rm{phy}}\leqslant 98.92\%$
in our simulations. In particular, models with high accretion rates yield
relatively low efficiencies of $64.9\%\leqslant\eta^{\rm{phy}}\leqslant
92.81\%$, because the hydrogen dissociation consumes a significant amount of
energy and this latent heat is retained in the post-shock region by H,
lowering $\eta^{\rm{phy}}$. In contrast, the low accretion rate models yield
that $91.85\%\leqslant\eta^{\rm{phy}}\leqslant 98.92\%,$ as molecular hydrogen
is largely retained in the pre-shock region, behaving like perfect gas.
Without surprise, $\eta^{\rm{phy}}$ v.s. $\mathcal{M}_{\rm{pre}}$ (the pre-
shock gas Mach number) of the low accretion models is largely consistent with
Figure 4 of Marleau et al. (2019).

Figure 2: Global shock efficiency (left) and post-shock entropy (right). They
are the plots of $\eta^{\rm{phy}}$ v.s. $\mathcal{M}_{\rm{pre}}$, and
$s_{\rm{ps}}^{l8}$ v.s. $M_{p}$ in Table 1. The red and black colors indicate
$10^{-2}$M${}_{\oplus}\cdot\rm{yr}^{-1}$ and
$10^{-3}$M${}_{\oplus}\cdot\rm{yr}^{-1}$accretion rate. The filled circles and
squares are the results of model 3, 5, 9, and 11 which use smaller
$r_{\rm{in}}$.

###  3.3 Post-shock entropy

Incorporating all simulation results, we present $s_{\rm ps}$ v.s. $M_{p}$ in
Figure 2. We find a clear dichotomy, characterized by a high-entropy “hot
group” with $s_{\rm ps}\gtrsim 16k_{\rm{B}}/m_{\ce{H}}$ (marked in orange),
and a low-entropy “cold group” with $s_{\rm ps}\lesssim
12k_{\rm{B}}/m_{\ce{H}}$ (marked in blue). By noting the fraction of H in the
pre-shock and post-shock regions in Table 1, we clearly identify that the “hot
group” is associated with models where most H2 is dissociated before reaching
the shock, whereas those without H2 dissociation before the shock fall into
the “cold group”. We should notice that the post-shock gas will continue to
cool in a radiative zone, further decreasing the entropy, and the general
zero-age entropy of planets (Spiegel & Burrows, 2012) should be lower than
$s_{\rm{ps}}$.

A few comments are in order concerning the dichotomy. First, the post-shock
entropy generally increases monotonically with increasing accretion rate and
planet mass. When choosing intermediate accretion rates for low mass planets
$M_{p}=$0.3-1MJ, one would fill the gap of $s_{\rm ps}$ in Figure 2 between
the two groups. This also applies to model 11 and 12 for massive planets with
low accretion rates. On the other hand, the dichotomy remains a valid
description, as the “gap” region corresponds to H2 being partially
dissociated. Second, Marleau et al. (2019) varied $\gamma$ from 1.1 to 1.667
to mimic realistic gas. Their entropy profiles with low $\gamma$ are similar
to our results when H2 dissociates and their high $\gamma$ results are similar
to our results when H2 does not dissociate. Third, even though $s_{\rm{ps}}$
of massive planets ($>3$MJ) can be low at low accretion rates
($<10^{-3}$M${}_{\oplus}\cdot\rm{yr}^{-1}$), we should still expect high
initial planet entropy because they must build the bulk of their mass through
phases of high accretion rates.

Finally, given its importance, we provide an approximate estimate of the
condition for which H2 gets fully dissociated. This can be obtained by
equating radiation temperature from the accretion luminosity at the pre-shock
region to $T_{\rm{g,dis}}$. During dissociation, we have
$F_{r}=cf_{r}E_{r}=cf_{r}a_{\rm{R}}T_{\rm{rad}}^{4}$ where $f_{r}\approx 1$.
Assuming the shock luminosity is much larger than the internal luminosity, we
obtain

$\displaystyle\eta^{\rm{phy}}\frac{GM_{p}\dot{M_{p}}}{r_{p}}$
$\displaystyle\approx$ $\displaystyle 4\pi
f_{c}r_{p}^{2}cf_{r}a_{\rm{R}}T_{\rm{g,dis}}^{4},$ (19) $\displaystyle\approx$
$\displaystyle 6\times
10^{-4}f_{c}f_{r}\bigg{(}\frac{r_{p}}{R_{J}}\bigg{)}^{2}L_{\odot}.$ (20)

We can see that the aforementioned $M_{p}/r_{p}$ is correlated to the
dissociation of H2. Equation 19 can be translated to a critical accretion rate
for H2 dissociation, given by

$\dot{M}_{\rm{crit}}\approx 7\times
10^{-4}\frac{f_{c}f_{r}}{\eta^{\rm{phy}}}\bigg{(}\frac{r_{p}}{R_{J}}\bigg{)}^{3}\bigg{(}\frac{M_{p}}{M_{J}}\bigg{)}^{-1}\frac{M_{\oplus}}{\rm
yr}.$ (21)

We can calculate $\dot{M}_{\rm{crit}}$ by plugging in our assumption of
$f_{c}=1$, and $f_{r}=1$ with $\eta^{\rm phy}$ from the simulation results. We
present $\dot{M}_{p}/\dot{M}_{\rm{crit}}$ in Table 1. When
$\dot{M}_{p}/\dot{M}_{\rm{crit}}>1$, H2 is indeed almost fully dissociated;
when $\dot{M}_{p}/\dot{M}_{\rm{crit}}\lessapprox 1$, H2 is partially
dissociated; when $\dot{M}_{p}/\dot{M}_{\rm{crit}}<0.56$, H2 largely remains
in molecular form in the pre-shock region. In general, we find that this
criterion well characterizes the chemical state of the hydrogen in the pre-
shock region and separates the “cold group” and “hot group” in our simulations
(see Figure 2).

##  4 Summary and discussion

In this letter, we carried out 12 radiation hydrodynamic simulations in 1D
with hydrogen EoS to study planetary accretion with a broad range of planet
masses and accretion rates. The simulations suggest that hydrogen dissociation
in the pre-shocked infalling gas plays an important role in setting the
initial condition of giant planets. More specifically, we find

  1. 1.

The post-shock entropy $s_{\rm ps}$ of the 12 simulations generally fall into
two groups (Figure 2): a “cold” group with low degree of H2 dissociation and a
“hot” group with high degree of H2 dissociation.

  2. 2.

The global shock efficiency $\eta^{\rm{phy}}$ \- the fraction of the accretion
energy that is radiated away - can be lowered by hydrogen dissociation.

  3. 3.

There is a critical accretion rate $\dot{M}_{\rm{crit}}$ above which H2 gets
largely dissociated in the infalling gas, which is given by Equation 21.

We note that the post-shock entropy (Figure 2) should be higher than the
entropy of zero-age planets because the newly accreted gas can still radiate
energy away. In addition, our simulations represent individual snapshots of
planetary accretion shocks, while the final result should depend on the
accretion history, especially over the period when giant planets build most of
their masses. Our results of the global shock efficiencies and post-shock
entropy can serve as more realistic inputs for more detailed modeling of the
structure and evolution of accreting planets (Mordasini et al., 2012; Berardo
et al., 2017), which eventually determines the initial conditions of gas
giants.

As a first approximation, we have assumed that accretion is spherically
symmetric. While we expect the actual accretion process could be more complex
through a CPD (Takasao et al., 2021), we note that existing studies of CPDs
have yet to self-consistently incorporate the major physical ingredients
including radiation with realistic EoS and magnetic fields. Despite the major
uncertainties in our ignorance of how gas giants accrete, our results also
imply that whether the hydrogen is accreted in molecular or atomic form likely
has significant impact in setting the initial condition of gas giants.

Current imaging surveys may be approaching the limit to detect more young gas
giants in PPDs, especially in systems with large inner cavities (Asensio-
Torres et al., 2021, but see Sanchis et al., 2020). There are also giant
planets with predicted dynamical mass awaiting for direct imaging detection
(e.g., Dong et al., 2015; Maire et al., 2017; Brown-Sevilla et al., 2021),
which are ideal tools to test the evolution models of gas giants. The
discovery space will be greatly enhanced with the upcoming James Webb Space
Telescope (JWST, Carter et al., 2021), the Extremely Large Telescope (ELT,
Carlomagno et al., 2020), and the Chinese Space Station Telescope (CSST)
which, together with modeling effort, will likely yield a more decisive
picture of gas giant formation.

We are indebted to Gabriel-Dominique Marleau for valuable comments on a
preliminary version of the paper, and thank Yuhiko Aoyama, Ruobing Dong, Chris
Ormel, and Shude Mao for helpful discussions. We also thank the anonymous
referee whose suggestions improved the quality of this Letter. This work is
supported by the National Key R&D Program of China (No.2019YFA0405100). ZC is
grateful to the CITA National Postdoctoral Fellowship, Tsinghua Astrophysics
Outstanding Fellowship, and Shuimu Tsinghua Scholar Program.

## References

  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41, doi: 10.3847/2041-8213/aaf741
  * Aoyama & Ikoma (2019) Aoyama, Y., & Ikoma, M. 2019, ApJ, 885, L29, doi: 10.3847/2041-8213/ab5062
  * Aoyama et al. (2018) Aoyama, Y., Ikoma, M., & Tanigawa, T. 2018, ApJ, 866, 84, doi: 10.3847/1538-4357/aadc11
  * Asensio-Torres et al. (2021) Asensio-Torres, R., Henning, T., Cantalloube, F., et al. 2021, A&A, 652, A101, doi: 10.1051/0004-6361/202140325
  * Avenhaus et al. (2018) Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44, doi: 10.3847/1538-4357/aab846
  * Ayliffe & Bate (2009) Ayliffe, B. A., & Bate, M. R. 2009, MNRAS, 397, 657, doi: 10.1111/j.1365-2966.2009.15002.x
  * Bae et al. (2021) Bae, J., Teague, R., & Zhu, Z. 2021, ApJ, 912, 56, doi: 10.3847/1538-4357/abe45e
  * Bae & Zhu (2018) Bae, J., & Zhu, Z. 2018, ApJ, 859, 118, doi: 10.3847/1538-4357/aabf8c
  * Balay et al. (1997) Balay, S., Gropp, W. D., McInnes, L. C., & Smith, B. F. 1997, in Modern Software Tools in Scientific Computing, ed. E. Arge, A. M. Bruaset, & H. P. Langtangen (Birkhäuser Press), 163–202
  * Balay et al. (2019) Balay, S., Abhyankar, S., Adams, M. F., et al. 2019, PETSc Users Manual, Tech. Rep. ANL-95/11 - Revision 3.11, Argonne National Laboratory
  * Berardo et al. (2017) Berardo, D., Cumming, A., & Marleau, G.-D. 2017, ApJ, 834, 149, doi: 10.3847/1538-4357/834/2/149
  * Brown-Sevilla et al. (2021) Brown-Sevilla, S. B., Keppler, M., Barraza-Alfaro, M., et al. 2021, A&A, 654, A35, doi: 10.1051/0004-6361/202140783
  * Carlomagno et al. (2020) Carlomagno, B., Delacroix, C., Absil, O., et al. 2020, Journal of Astronomical Telescopes, Instruments, and Systems, 6, 035005, doi: 10.1117/1.JATIS.6.3.035005
  * Carter et al. (2021) Carter, A. L., Hinkley, S., Bonavita, M., et al. 2021, MNRAS, 501, 1999, doi: 10.1093/mnras/staa3579
  * Chen & Bai (in prep) Chen, Z., & Bai, X. in prep
  * Chen et al. (2019) Chen, Z., Coleman, M. S. B., Blackman, E. G., & Frank, A. 2019, Journal of Computational Physics, 388, 490, doi: 10.1016/j.jcp.2019.03.016
  * de Val-Borro et al. (2007) de Val-Borro, M., Artymowicz, P., D’Angelo, G., & Peplinski, A. 2007, A&A, 471, 1043, doi: 10.1051/0004-6361:20077169
  * Dipierro et al. (2015) Dipierro, G., Price, D., Laibe, G., et al. 2015, MNRAS, 453, L73, doi: 10.1093/mnrasl/slv105
  * Dong et al. (2017) Dong, R., Li, S., Chiang, E., & Li, H. 2017, ApJ, 843, 127, doi: 10.3847/1538-4357/aa72f2
  * Dong et al. (2015) Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015, ApJ, 809, L5, doi: 10.1088/2041-8205/809/1/L5
  * Fung et al. (2019) Fung, J., Zhu, Z., & Chiang, E. 2019, ApJ, 887, 152, doi: 10.3847/1538-4357/ab53da
  * Goldreich & Tremaine (1979) Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857, doi: 10.1086/157448
  * Goodman & Rafikov (2001) Goodman, J., & Rafikov, R. R. 2001, ApJ, 552, 793, doi: 10.1086/320572
  * Gressel et al. (2013) Gressel, O., Nelson, R. P., Turner, N. J., & Ziegler, U. 2013, ApJ, 779, 59, doi: 10.1088/0004-637X/779/1/59
  * Haffert et al. (2019) Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nature Astronomy, 3, 749, doi: 10.1038/s41550-019-0780-5
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Keppler et al. (2018) Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44, doi: 10.1051/0004-6361/201832957
  * Levermore & Pomraning (1981) Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321, doi: 10.1086/159157
  * Lin & Papaloizou (1986) Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846, doi: 10.1086/164653
  * Long et al. (2018) Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17, doi: 10.3847/1538-4357/aae8e1
  * Maire et al. (2017) Maire, A. L., Stolker, T., Messina, S., et al. 2017, A&A, 601, A134, doi: 10.1051/0004-6361/201629896
  * Malygin et al. (2014) Malygin, M. G., Kuiper, R., Klahr, H., Dullemond, C. P., & Henning, T. 2014, A&A, 568, A91, doi: 10.1051/0004-6361/201423768
  * Marleau et al. (2017) Marleau, G.-D., Klahr, H., Kuiper, R., & Mordasini, C. 2017, ApJ, 836, 221, doi: 10.3847/1538-4357/836/2/221
  * Marleau et al. (2019) Marleau, G.-D., Mordasini, C., & Kuiper, R. 2019, ApJ, 881, 144, doi: 10.3847/1538-4357/ab245b
  * Marley et al. (2007) Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007, ApJ, 655, 541, doi: 10.1086/509759
  * Mordasini et al. (2012) Mordasini, C., Alibert, Y., Georgy, C., et al. 2012, A&A, 547, A112, doi: 10.1051/0004-6361/201118464
  * Omukai et al. (2005) Omukai, K., Tsuribe, T., Schneider, R., & Ferrara, A. 2005, ApJ, 626, 627, doi: 10.1086/429955
  * Pérez et al. (2016) Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Science, 353, 1519, doi: 10.1126/science.aaf8296
  * Pinte et al. (2018) Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJ, 860, L13, doi: 10.3847/2041-8213/aac6dc
  * Pinte et al. (2020) —. 2020, ApJ, 890, L9, doi: 10.3847/2041-8213/ab6dda
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62, doi: 10.1006/icar.1996.0190
  * Sanchis et al. (2020) Sanchis, E., Picogna, G., Ercolano, B., Testi, L., & Rosotti, G. 2020, MNRAS, 492, 3440, doi: 10.1093/mnras/staa074
  * Semenov et al. (2003) Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611, doi: 10.1051/0004-6361:20031279
  * Spiegel & Burrows (2012) Spiegel, D. S., & Burrows, A. 2012, ApJ, 745, 174, doi: 10.1088/0004-637X/745/2/174
  * Szulágyi & Ercolano (2020) Szulágyi, J., & Ercolano, B. 2020, ApJ, 902, 126, doi: 10.3847/1538-4357/abb5a2
  * Szulágyi et al. (2014) Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65, doi: 10.1088/0004-637X/782/2/65
  * Szulágyi & Mordasini (2017) Szulágyi, J., & Mordasini, C. 2017, MNRAS, 465, L64, doi: 10.1093/mnrasl/slw212
  * Takasao et al. (2021) Takasao, S., Aoyama, Y., & Ikoma, M. 2021, ApJ, 921, 10, doi: 10.3847/1538-4357/ac0f7e
  * Tanigawa et al. (2012) Tanigawa, T., Ohtsuki, K., & Machida, M. N. 2012, ApJ, 747, 47, doi: 10.1088/0004-637X/747/1/47
  * Teague et al. (2019) Teague, R., Bae, J., & Bergin, E. A. 2019, Nature, 574, 378, doi: 10.1038/s41586-019-1642-0
  * Thanathibodee et al. (2019) Thanathibodee, T., Calvet, N., Bae, J., Muzerolle, J., & Hernández, R. F. 2019, ApJ, 885, 94, doi: 10.3847/1538-4357/ab44c1
  * van der Marel et al. (2013) van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2013, Science, 340, 1199, doi: 10.1126/science.1236770
  * Wagner et al. (2018) Wagner, K., Follete, K. B., Close, L. M., et al. 2018, ApJ, 863, L8, doi: 10.3847/2041-8213/aad695
  * Zhu & Stone (2014) Zhu, Z., & Stone, J. M. 2014, ApJ, 795, 53, doi: 10.1088/0004-637X/795/1/53

##  Appendix A Equation of state

The general EoS can be derived based on the partition functions of the
underlying species. The partition functions of all the species used in this
work are given by

$\displaystyle Z_{\ce{H2}}$ $\displaystyle=$ $\displaystyle\left(\frac{2\pi
m_{\ce{H2}}k_{\rm{B}}T_{\rm{g}}}{h^{2}}\right)^{3/2},$ (A1) $\displaystyle
Z_{\ce{H}}$ $\displaystyle=$ $\displaystyle\left(\frac{2\pi
m_{\ce{H}}k_{\rm{B}}T_{\rm{g}}}{h^{2}}\right)^{3/2}\exp\left(-\frac{\phi_{\rm{dis}}}{2k_{\rm{B}}T_{\rm{g}}}\right),$
(A2) $\displaystyle Z_{\ce{H+}}$ $\displaystyle=$
$\displaystyle\left(\frac{2\pi
m_{\ce{H+}}k_{\rm{B}}T_{\rm{g}}}{h^{2}}\right)^{3/2}\exp\left(-\frac{\phi_{\rm{dis}}+2\phi_{\rm{ion}}}{2k_{\rm{B}}T_{\rm{g}}}\right),$
(A3) $\displaystyle Z_{\ce{e-}}$ $\displaystyle=$
$\displaystyle\left(\frac{2\pi
m_{\ce{e-}}k_{\rm{B}}T_{\rm{g}}}{h^{2}}\right)^{3/2},$ (A4) $\displaystyle
Z_{\ce{He}}$ $\displaystyle=$ $\displaystyle\left(\frac{2\pi
m_{\ce{He}}k_{\rm{B}}T_{\rm{g}}}{h^{2}}\right)^{3/2},$ (A5)

where $m_{\ce{H2}}$, $m_{\ce{H}}$, $m_{\ce{H+}}$, $m_{\ce{e-}}$, and
$m_{\ce{He}}$ are the masses of the elements. $\phi_{\rm{dis}}=7.17\times
10^{-12}$ erg and $\phi_{\rm{ion}}=2.18\times 10^{-11}$ erg are the binding
energy of H2 and H.

In the above, we do not consider the spin and atomic structure of any species,
and have ignored the rotational and vibrational degrees of freedom of H2.
Incorporating these degrees of freedom would lead to much complex partition
functions, yet because the energy needed to excite these transitions are small
compared to the binding energy, we expect the resulting impact to be
relatively minor. Moreover, we need to use a tabulated EoS if we were to
include them to the partition function, which would make our calculations
inefficient and less accurate.

In many of the simulations, H2 reforms in the post-shock region. We note that
the underlying assumption with our treatment of the general EoS is chemical
equilibrium, implying that H2 formation is instantaneous based on the Saha
equation. In reality, given that the post-shock gas density is high
($\rho_{\rm{ps}}>10^{-8}$ g$\cdot$cm-3), we anticipate that efficient H2
formation can be achieved through the three-body H2 formation channel on a
timescale of less than one minute (Omukai et al., 2005), which is much shorter
than the free-fall time.

##  Appendix B Global shock efficiency

Global shock efficiency $\eta^{\rm{phy}}$ quantifies the fraction of the
energy that leave the system, in our case the radiation energy. For steady
state solutions, Equation 3 and 4 become,

$\displaystyle\frac{1}{r^{2}}\partialderivative{r}[r^{2}(\epsilon_{\rm
kin}+\epsilon+p)v]$ $\displaystyle=$ $\displaystyle\mathbb{G}-\rho
v\nabla\phi,$ (B1)
$\displaystyle\frac{1}{r^{2}}\partialderivative{r}(r^{2}F_{r})$
$\displaystyle=$ $\displaystyle-\mathbb{G},$ (B2)

where $\phi=-GM_{p}/r$ is the gravitational potential energy. We can multiply
by $r^{2}$ and add up the two equations above, recognizing that
$\dot{M}_{p}=-4\pi\rho r^{2}v$ is constant in steady state and $4\pi
r^{2}F_{r}=L_{r}$ is the luminosity, to obtain

$\dot{M}_{p}\partialderivative{r}(\epsilon+\epsilon_{\rm{kin}}+p/\rho+\phi)+\partialderivative{L_{r}}{r}=0\
.$ (B3)

Equation B3 illustrates the balance of the radiative flux and the total energy
flux $\dot{E}(r)=-\dot{M}_{p}(\epsilon+\epsilon_{\rm{kin}}+p/\rho+\phi)$ in
steady state. Since radiation energy is the one that leaves the system, the
fraction of accretion energy radiated away, which is defined as the shock
efficiency $\eta^{\rm{phy}}$, is given by

$\eta^{\rm{phy}}=\frac{L_{\rm{r}}(r_{\rm{max}})-L_{\rm{r}}(r_{\rm{ps}})}{\dot{E}(r_{\rm{max}})}=\frac{\dot{E}(r_{\rm{max}})-\dot{E}(r_{\rm{ps}})}{\dot{E}(r_{\rm{max}})}\leq
1\ .$ (B4)

Consequently, $1-\eta^{\rm{phy}}$ is the fraction of the energy retained after
the shock.


# Work distributions in the $T=0$ Random Field Ising Model

Xavier Illa  Department of Applied Physics, Helsinki University of Technology,
P.O.Box 1100, FIN-02015 HUT, Finland  Josep Maria Huguet  Departament de
Física Fonamental, Universitat de Barcelona
Martí i Franquès 1, Facultat de Física, 08028 Barcelona, Catalonia  Eduard
Vives  Departament d’Estructura i Constituents de la Matèria, Universitat de
Barcelona
Martí i Franquès 1 , Facultat de Física, 08028 Barcelona, Catalonia

###### Abstract

We perform a numerical study of the three-dimensional Random Field Ising Model
at $T=0$. We compare work distributions along metastable trajectories obtained
with the single-spin flip dynamics with the distribution of the internal
energy change along equilibrium trajectories. The goal is to investigate the
possibility of extending the Crooks fluctuation theorem G.E.Crooks (1999) to
zero temperature when, instead of the standard ensemble statistics, one
considers the ensemble generated by the quenched disorder. We show that a
simple extension of Crooks fails close to the disordered induced equilibrium
phase transition due to the fact that work and internal energy distributions
are very asymmetric.

###### pacs:

75.60.Ej, 75.50.Lk, 81.30.Kf, 81.40.Jj

††preprint: EV/25-07-2008

##  I Introduction

The $T=0$ Random Field Ising Model is a prototype model for the study of
collective phenomena in disordered systems. Although it neglects thermal
fluctuations, it contains essential competitions between the quenched
disorder, the ferromagnetic interaction and the external applied field. The
model can be numerically studied from two different points of view: on the one
hand, the exact ground state calculation Hartmann and Young (2001); Middleton
and D.S.Fisher (2002); Dukovski and Machta (2003); Wu and Machta (2005)
provides an approach to the equilibrium phase diagram. On the other the use of
a local relaxation dynamics based on single spin-flips provides a good
framework for the understanding of avalanche dynamics and hysteresis Sethna et
al. (1993); J.P.Sethna et al. (2005), which is closer to experimental
observations. In this sense, the model is a good workbench for the comparison
of equilibrium and out-of-equilibrium trajectories.

A number of non-equilibrium work theorems Sevick et al. (2008) have received a
lot of attention in the last 10 years, particularly after the work of
Jarzynski in 1997 C.Jarzynski (1997); Jarzynski (2006). These theorems relate
in different ways the distribution of work performed on a system which is
driven out of equilibrium to some equilibrium thermodynamic properties. One
example is the original Jarzinsky’s equality: $\langle e^{-\beta
W}\rangle=e^{-\beta\Delta F}$, where $W$ is the work performed on the system
that has been driven (out-of-equilibrium) by varying an external control
parameter $H$ changing from $H(0)$ to $H(1)$, and $\Delta F$ is the free
energy difference between two states $0$ and $1$ that correspond to the
equilibrium states at $H(0)$ and $H(1)$. The average $\langle\cdot\rangle$
should be understood as obtained after many repetitions of the driving
process. The system is assumed to be in contact with a heat bath at
temperature $\beta^{-1}$ which generates an statistical ensemble of copies of
the system which is the source of work fluctuations. Most of these theorems
have been experimentally verified. Collin et al. (2005); Trepagnier et al.
(2004)

The goal of this paper is to investigate how such kind of theorems can be
extended to systems at $T=0$. In such a case, although equilibrium states and
out-of-equilibrium trajectories are well defined, there are no thermal
fluctuations. A priori it seems that there is no statistical ensemble over
which one can define probability distributions or averages. The idea we want
to test is whether or not, for systems with quenched disorder, the thermal
ensemble can be substituted by the ensemble of different realizations of
disorder and still some work theorems can be applied.

In order to perform the $T\rightarrow 0$ limit we choose a different but
related work theorem which is Crooks fluctuations theorem G.E.Crooks (1999);
Crooks (2000). The advantage will be that it allows to derive an equality that
can be extrapolated to the $T\rightarrow 0$ limit. This theorem can be written
as:

$\frac{P_{F}(W)}{P_{R}(-W)}=e^{\beta(W-\Delta F)}$ (1)

In this case $P_{F}(W)$ and $P_{R}(-W)$ are the probability densities
corresponding to the out-of-equilibrium work performed on a system driven
forward from $H(0)$ to $H(1)$ and reversely from $H(1)$ to $H(0)$. The Crooks
fluctuation theorem has been formulated under several assumptions: the driven
systems must be finite, classical and coupled to a thermal bath. Moreover, the
dynamics should be stochastic, Markovian and reversible and the entropy
production should be odd under time reversal. Some of these assumptions are
clearly not accomplished at the work at hand. For instance the thermal bath is
at $T=0$ and the dynamics is deterministic. Nevertheless we have been a little
speculative and investigated the possibility of extending the Crooks
fluctuation theorem to $T=0$ for systems with quenched disorder.

From Eq. (1) one can derive that the value $W^{*}$ for which the two
probability densities are equal, i.e. $P_{F}(W^{*})=P_{R}(-W^{*})$ satisfies
$W^{*}=\Delta F$. This result is particularly suitable to be investigated at
$T\rightarrow 0$. Note that in the low temperature limit, the only possibility
to have a crossing point of the probability densities at $W^{*}$ is that
$W^{*}=\Delta F(T\rightarrow 0)=\Delta U$, so that the diverging behaviour of
$\beta$ might be cancelled to get a finite limit
$P_{F}(W^{*})/P_{R}(-W^{*})\rightarrow 1$.

##  II Model

For our investigation we consider the Random Field Ising Model. A set of $N$
spin variables $S_{i}=\pm 1$ with $i=1,\cdots N$ is defined on a regular 3D
cubic lattice with linear size $L$ ($N=L\times L\times L$). We are interested
in following the response of the order parameter (magnetization)
$M=\sum_{i=1}^{N}S_{i}$ as a function of the external driving field $H$. The
Hamiltonian (magnetic enthalpy) of the model is defined as:

${\cal H}=U-HM$ (2) $U=-\sum_{ij}S_{i}S_{j}-\sum_{i=1}^{N}h_{i}S_{i}$ (3)

where the first sum in Eq. 3, extending over nearest neigbours pairs, accounts
for a ferromagnetic exchange interaction, the second term in Eq. 3 includes
the interaction with the quenched disorder and the second term in Eq. 2
accounts for the interaction with the external field $H$ that will be used as
the driving parameter. The local fields $h_{i}$ are independent random
variables, quenched, Gaussian distributed with zero mean and standard
deviation $\sigma$. All along the paper we will use small letters
($u=\frac{U}{N}$,$m=\frac{M}{N}\ldots$) to refer to intensive magnitudes.

The equilibrium properties of this model can be studied at $T=0$. For a given
realization of the random fields $\\{h_{i}\\}$ and a given value of $H$, the
exact ground state can be found in a polynomial time Middleton (2001).
Moreover optimized algorithms can also be used to obtain the full sequence of
ground-sates as a function of $H$ between the two saturated states
$\\{S_{i}=-1\\}$ at $H=-\infty$ and $\\{S_{i}=1\\}$ at $H=\infty$ C.Frontera
et al. (2000). In the thermodynamic limit the 3D model exhibits a first order
phase transition at $H=0$ and for $\sigma<\sigma_{c}^{eq}=2.27$ Hartmann and
Young (2001); Middleton and D.S.Fisher (2002). Although finite systems do not
exhibit a true phase transition, there is a region of $\sigma$ and $H$ where
the correlation length is of the same order as the system size. This leads to
collective effects and correlations involving all the spins of the system.

Concerning the non-equilibrium trajectories, the model has been intensively
studied by using the so called single-spin flip metastable dynamics at $T=0$
which can be understood as a $T\rightarrow 0$ Glauber dynamics. It consists in
adiabatically sweeping the external field and relaxing individual spins
according to their local energy, i.e. spins should align with the local field

$F_{i}=H+h_{i}+\sum_{k}S_{k}$ (4)

where the sum extends over the neigbours of spin $i$. Note that the reversal
of a spin at a given field $H$ may induce the reversal of some of its
neigbours at the same field $H$. Such collective events are the so called
avalanches that end when all the spins are stable. Only after the avalanches
have finished, the external field is varied again. It has been shown that,
during a monotonous driving of the external field, no reverse spin-flips may
occur and, moreover, the model exhibits the so called abelian property: i.e.
the final states do not depend on the order in which the unstable spins are
flipped. In the thermodynamic limit, the system with this metastable dynamics
also shows a critical point at $\sigma_{c}^{met}=2.21$ Sethna et al. (1993);
F.J.Pérez-Reche and E.Vives (2003), below which the metastable trajectory
exhibits a discontinuity.

Note that avalanches are the source of energy dissipation Ortín and Goicoechea
(1998); X.Illa et al. (2005). If we think about an increasing field
trajectory, the triggering spin of each avalanche flips at the field $H$ that
corresponds to $F_{i}=0$. This spin, therefore, flips without energy loss
$\Delta{\cal H}=\Delta U-H\Delta M=0$. But the subsequent unstable spins flip
(at a fixed external field $H$) with an external force $F_{i}>0$ giving raise
to an energy loss (associated to each individual spin flip) $Q=-2F_{i}<0$ with
$Q=\Delta U-W$. For the decreasing field branches, $F_{i}<0$ but the energy
loss is then $Q=2F_{i}<0$ since the spins flip from $1$ to $-1$.

Note that in this discussion we are using the standard definition of work that
is not a state function. It is computed as a sum over the ($k=1,2\ldots$)
spins that flip along the trajectory

$W=\sum_{k}H_{k}\Delta M_{k},$ (5)

where $\Delta M_{k}=2$. Recently there has been a discussion Vilar and Rubi
(2008); A.Imparato and L.Peliti (2007); J.Horowitz and C.Jarzynski (2008);
Peliti (2008) on which is the most suitable definition of work to be used in
such non-equilibrium work theorems. Already from the initial Jarzynski works
C.Jarzynski (1997), it was proposed that if the system is driven by
controlling the external field, the convenient definition of work is the
integral over the trajectory $V=\int MdH$. Note that along a metastable
trajectory , the two definitions of work are related by $W=\Delta(HM)-V$.
Without going into the discussion, we will numerically test here the two
possibilities. In Sec. III we will test whether the crossing point $W^{*}$ of
the histograms $P_{F}(W)$ and $P_{R}(-W)$ satisfies

$W^{*}=\langle\Delta U\rangle$ (6)

and, in Sec. IV we will test whether the crossing point $V^{*}$ of the
histograms $P_{F}(V)$ and $P_{R}(-V)$ satisfies the corresponding equation:

$V^{*}=\langle\Delta U-\Delta(HM)\rangle$ (7)

The test of the two hypothesis will require slightly different strategies (as
indicated in Fig. 1) since $V$ can not be computed for trajectories starting
at saturation ($H\rightarrow\pm\infty$).

Figure 1:  Schematic representation of the two different strategies that have
been used. (a) corresponds to the first strategy in which the starting point
is the saturated state $0$ ($m=-1$,$H=-\infty$), and (b) corresponds to the
second strategy, starting from the equilibrium state $0$ at a finite field
$H_{0}$. Dashed lines correspond to metastable trajectories and the continuum
line indicates the equilibrium states.

In the first strategy (Sec. III), we will compute forward trajectories
starting from the state $0$, saturated with $m=-1$ at $H=-\infty$ and sweep
the field adiabatically from $-\infty$ to $H_{1}$ until the metastable state
$A$ is reached. Then we will compute the equilibrium state $1$ corresponding
to the field $H_{1}$ and perform the reverse trajectory from $1$ to $0$ with
the metastable dynamics. We will compute $\Delta u=(U_{1}-U_{0})/N$ and the
following works (using the ’standard’ definition):

$\displaystyle w_{0\rightarrow 1}$ $\displaystyle=$
$\displaystyle\frac{W_{0\rightarrow
1}}{N}=\frac{1}{N}\int_{0}^{A}HdM+\frac{1}{N}\int_{A}^{1}HdM=\frac{1}{N}\int_{0}^{A}HdM+\frac{1}{N}H(M_{1}-M_{A})$
(8) $\displaystyle w_{1\rightarrow 0}$ $\displaystyle=$
$\displaystyle\frac{W_{1\rightarrow 0}}{N}=\frac{1}{N}\int_{1}^{0}HdM$ (9)

where $M_{A}$ and $M_{1}$ are the magnetizations of states $A$ and $1$. The
integrals are computed along the trajectories schematically represented in
Fig. 1(a) using Eq. 5. Note also that the second integral in Eq. 8 can be
computed since the process $A\rightarrow 1$ takes place at constant field
$H_{1}$.

In the second strategy (Sec. IV), we will start from a computed equilibrium
state at $H_{0}$, perform a metastable trajectory until the state $A$ is
reached at $H_{1}$. Then we will compute the equilibrium state $1$ at $H_{1}$
and perform a metastable trajectory driving back until reaching the state $B$
at $H_{0}$. In this case we will use the alternative definition of work:

$\displaystyle v_{0\rightarrow 1}=v_{0\rightarrow A}=\frac{V_{0\rightarrow
1}}{N}$ $\displaystyle=$ $\displaystyle\frac{V_{0\rightarrow
A}}{N}=\frac{1}{N}\int_{0}^{A}MdH,$ (10) $\displaystyle v_{1\rightarrow
0}=v_{1\rightarrow B}=\frac{V_{1\rightarrow 0}}{N}$ $\displaystyle=$
$\displaystyle\frac{V_{1\rightarrow B}}{N}=\frac{1}{N}\int_{1}^{B}MdH.$ (11)

The integrals are computed along the trajectories schematically represented in
Fig. 1(b) using $V=\sum_{k}M_{k}\Delta H_{k}$ where $k$ accounts for the
sequence of interavalanche field increments $\Delta H_{k}$ occuring at
constant $M_{k}$ along the trajectory.

.

Figure 2:  (Color on line) Example of histograms corresponding to
$P_{F}(w_{0\rightarrow 1})$, $P_{R}(-w_{1\rightarrow 0})$ and $P(\Delta u)$.
The arrow pointing downwards indicates the position of the maximum of
$P(\Delta u)$ and the arrow pointing upwards the average value $\langle\Delta
u\rangle$. The example corresponds to a system with $\sigma=4$, $H_{1}=-0.35$
and $L=12$. Histograms have been obtained by cumulating data corresponding to
$4\times 10^{5}$ realizations of disorder.

##  III First strategy

Fig. 2 shows histograms corresponding to the distributions of $w_{0\rightarrow
1}$, $\Delta u$ and $-w_{1\rightarrow 0}$ for $\sigma=4$, $H_{1}=-0.35$ and a
system size $L=12$. The computed histograms give an estimation of the
corresponding probability densities $P_{F}(w_{0\rightarrow 1})$,
$P_{R}(-w_{1\rightarrow 0})$ and $P(\Delta u)$. Four important points should
be realized from this example: (i) For this value of $\sigma$ and field
$H_{1}$, the densities look symmetric and Gaussian. This will not be the case
when $H_{1}$ and $\sigma$ are close to the disorder induced phase transition
at $\sigma_{c}^{eq}$ and $H_{1}=0$. (ii) Second, the distribution of $\Delta
u$ has a width similar to those of the out-of-equilibrium works. This is
different from what happens in the analysis of Crooks fluctuation theorem at
finite $T$. Typically, equilibrium thermal fluctuations are much smaller than
work fluctuations. (iii) For increasing systems size the histograms become
narrower and then the crossing points are harder to locate (iv) Finally, note
that for this case the hypothesis that we are testing is fulfilled: the peak
in $P(\Delta u)$ as well as the average $\langle\Delta u\rangle$ coincide with
the crossing point $w^{*}$ of $P_{F}(w_{0\rightarrow 1})$ and
$P_{R}(-w_{1\rightarrow 0})$.

Figure 3:  (Color on line) Example of histograms corresponding to
$P_{F}(w_{0\rightarrow 1})$, $P_{R}(-w_{1\rightarrow 0})$ and $P(\Delta u)$.
The arrow pointing downwards indicates the position of the maximum of
$P(\Delta u)$ and the arrow pointing upwards the average value $\langle\Delta
u\rangle$. The example corresponds to a system with $\sigma=3$, $H_{1}=0.1$
and $L=12$. The histograms have been computed by cumulating data corresponding
to $1.6\times 10^{6}$ realizations of disorder. The inset shows a detailed
view of the crossing point of the $P_{F}(w_{0\rightarrow 1})$ and
$P_{R}(-w_{1\rightarrow 0})$ histograms. Figure 4:  Schematic representation
of the case $\sigma\lesssim\sigma_{c}^{eq}$ and $H_{1}>0$. Dashed lines
correspond to metastable trajectories, continuous lines to the equilibrium
trajectory and the dashed-dotted line to the non typical equilibrium states
that are responsible for the high asymmetry of the distribution
$P_{R}(-w_{1\rightarrow 0})$.

Fig. 3 shows histograms corresponding to $P(w_{0\rightarrow 1})$, $P(\Delta
u)$ and $P(-w_{1\rightarrow 0})$ for $\sigma=3$, $H_{1}=0.1$ and a system size
$L=12$. Note than in this case, the work distributions as well as the energy
distribution are very asymmetric. This is because although in this case we
expect $M_{1}>0$ since $H_{1}>0$, a certain non vanishing fraction of the
equilibrium states still has $M_{1}<0$, as schematically indicated in Fig.4.
In other words, for certain realizations of disorder, the state $1$ turns out
to be non-typical and previous to the equilibrium transition from $M<0$ to
$M>0$. Consequently, the distribution of reverse works $w_{1\rightarrow 0}$
widens enormously. As can be seen the proposed equality $w^{*}=\Delta u$
clearly fails in this case. The relative error is $|w^{*}-\langle\Delta
u\rangle|/|\langle\Delta u\rangle|\sim 0.25$.

Figure 5:  (Color on line) Comparison of the crossing point of the histograms
$w^{*}$ ($\Box$) with $\langle\Delta u\rangle$ ($\Circle$) as a function of
$H_{1}$, corresponding to (a) $\sigma=4$ and (b) $\sigma=3$ for different
values of $L$ as indicated by the legend. For the case $L=12$ and $\sigma=3$
the dashed line indicated the position of the maximum of $P(\Delta u)$. Data
for $L=18$ ($L=24$) has been shifted 0.3 (0.6) units along the vertical axis
in order to clarify the picture.

We have performed a detailed study of the deviation of $w^{*}$ from
$\langle\Delta u\rangle$ for different values of $\sigma$, $H_{1}$ and $L$.
Examples of the results are presented in Fig. 5. Note that the data
corresponding to the crossing point $w^{*}$ exhibits increasing error bars
with increasing $H_{1}$. This is because the histograms of
$P_{F}(w_{0\rightarrow 1})$, $P_{R}(-w_{1\rightarrow 0})$ become more and more
separated and the finding of the intersection requires more and more
statistics. (This problem is accentuated for larger values of $H_{1}$ and for
larger system sizes). One can conclude therefore than the proposed
extrapolation of Crooks theorem, given by the equality $w^{*}=\langle\Delta
u\rangle$ fails in the region of $\sigma$ and $H_{1}$ where the system
exhibits a collective behaviour (i.e. the correlation length is similar to the
system size) because of the proximity of the equilibrium phase transition.
Consistently, we observe in Fig. 5(b) (comparing the data corresponding to
$L=12$ and $L=18$ at $H_{1}\sim 0$) that the region of breakdown becomes
smaller when the system size $L$ is larger. As expected, increasing the system
size with fixed correlation length decreases the collective effects.

##  IV Second strategy

Fig. 6 shows histograms corresponding to the distributions of $v_{0\rightarrow
1}$ and $-v_{1\rightarrow 0}$ for the case $L=12$, $\sigma=3$, $H_{1}=-10$ and
$H_{2}=-0.4$. As can be seen the intersection point $v^{*}$ agrees very well
with the average value of $\Delta u-\Delta(Hm)$. In this case we have also
studied if this agreement is equally good in other regions of the phase
diagram.

Figure 6:  (Color on line) Example of histograms corresponding to
$P_{F}(v_{0\rightarrow 1})$, $P_{R}(-v_{1\rightarrow 0})$ and $P\left(\Delta
u-\Delta(Hm)\right)$. Data corresponds to $H_{0}=-10.0$,$H_{1}=-0.4$,
$\sigma=3$ and $L=12$. Histograms are computed by cumulating $10^{5}$
realizations of disorder. The arrows pointing down (up) indicate the peak
(average value) of the distribution $P(\Delta u-\Delta(HM))$.  Figure 7:
(Color on line) Comparison of the crossing point $v^{*}$ ($\Box$) of the
histograms $P_{F}(v_{0\rightarrow 1})$ and $P_{R}(-v_{1\rightarrow 0})$ with
$\langle\Delta u-\Delta(Hm)\rangle$ ($\Circle$) as a function of $H_{2}$,
corresponding to $\sigma=3$ and $L=12$. Data corresponds to averages over
$10^{5}$ realizations of disorder

Fig 7 shows a comparison of $\langle\Delta u\rangle$ and
$v^{*}+\langle\Delta(Hm)\rangle$. On the right hand side, the value of $v^{*}$
has not been obtained by a direct location of the intersection of the
histograms but locating the mid point between the positions of the maxima of
the histograms corresponding to $P_{F}(v_{0\rightarrow 1})$ and
$P_{R}(-v_{1\rightarrow 0})$ which are far separated but very symmetric. As
can be seen, as occurs with strategy 1, the agreement
$v^{*}\simeq\langle\Delta u-\Delta(Hm)\rangle$ fails when $H_{2}$ approaches
the phase transition region.

##  V Summary and conclusion

We have investigated the possibility of extending some work fluctuation
theorems to $T=0$ for systems with quenched disorder. In this case, thermal
averages should be substituted by disorder averages. We have proposed an
hypothesis based on the extrapolation of Crooks theorem to the $T\rightarrow
0$ limit and we have numerically tested its validity. The hypothesis is that
if one considers the distributions of non-equilibrium work corresponding to a
forward and backwards trajectory, the crossing point $w^{*}$ of
$P_{F}(w_{0\rightarrow 1})$ and $P_{R}(-w_{1\rightarrow 0})$, is equal to the
average of the internal energy difference $\Delta u=u_{1}-u_{2}$. The
investigation has been done by considering two strategies, based on the two
possible definitions of work that have been discussed in the literature:
$w=\int Hdm$ and $v=\int mdH$.

The reported numerical investigation indicates that, for both strategies, the
formulated hypothesis is valid when the system does not behave collectively.
Thus, far from the equilibrium critical point, where correlations are small
compared to system size, the distribution $p(\Delta u)$ is very symmetric and
the average (equilibrium) value $\langle\Delta u\rangle$ can be obtained from
the crossing point $w^{*}$ of the distributions of the non equilibrium works.

When the system is close to the critical point a definitive conclusion can not
be stated. Data for small $L$ suggest that when the correlation length
approaches the system size $L$, the distribution $p(\Delta u)$ is very
assymetric and there is no connection between the equilibrium value
$\langle\Delta u\rangle$ and the crossing point $w^{*}$ of the non-equilibrium
work distributions. A bigger computational effort (not affordable at present)
would be needed in order to compute the accurate histograms for $L>24$. The
presented results may provide some clues for future investigations in order to
extend work fluctuation theorems.

###### Acknowledgements.

This work has received financial support from CICyT (Spain), project
MAT2007-61200, CIRIT (Catalonia), project 2005SGR00969, EU Marie Curie RTN
MULTIMAT, Contract No. MRTN-CT-2004-5052226 and the Academy of Finland. X.I.
acknowledges the hospitality of ECM Department (Universitat de Barcelona)
where the work at hand was to a large degree completed. We also acknowledge
fruitful discussions with F. Ritort, A. Planes, and M.J. Alava.

## References

  * G.E.Crooks (1999) G.E.Crooks, Phys. Rev. E 60, 2721 (1999).
  * Hartmann and Young (2001) A. K. Hartmann and A. P. Young, Phys. Rev. B 64, 214419 (2001).
  * Middleton and D.S.Fisher (2002) A. A. Middleton and D.S.Fisher, Phys. Rev. B 65, 134411 (2002).
  * Dukovski and Machta (2003) I. Dukovski and J. Machta, Phys. Rev. B 67, 014413 (2003).
  * Wu and Machta (2005) Y. Wu and J. Machta, Physical Review Letters 95, 137208 (2005).
  * Sethna et al. (1993) J. P. Sethna, K. Dahmen, S. Kartha, J. A. Krumhansl, B. W. Roberts, and J. D. Shore, Phys. Rev. Lett. 70, 3347 (1993).
  * J.P.Sethna et al. (2005) J.P.Sethna, K.Dahmen, and O.Perković, _The Science of Hysteresis_ (Elsevier Inc., 2005), vol. 2, chap. 2, pp. 107–179.
  * Sevick et al. (2008) E. Sevick, R. Prabhakar, Stephen, and R. Williams, Annual Review of Physical Chemistry 59, 603 (2008).
  * C.Jarzynski (1997) C.Jarzynski, Phys. Rev. Lett. 78, 2690 (1997).
  * Jarzynski (2006) C. Jarzynski, Phys. Rev. E 73, 046105 (2006).
  * Collin et al. (2005) D. Collin, F. Ritort, C. Jarzynski, S. Smith, I. T. Jr, and C. Bustamante, Nature 437, 231 (2005).
  * Trepagnier et al. (2004) E. Trepagnier, C. Jarzynski, F. Ritort, G. Crooks, C. Bustamante, and J. Liphardt, PNAS 101, 15038 (2004).
  * Crooks (2000) G. E. Crooks, Phys. Rev. E 61, 2361 (2000).
  * Middleton (2001) A. A. Middleton, Phys. Rev. Lett. 88, 017202 (2001).
  * C.Frontera et al. (2000) C.Frontera, J.Goicoechea, J.Ortín, and E.Vives, J. Comp. Phys. 160, 117 (2000).
  * F.J.Pérez-Reche and E.Vives (2003) F.J.Pérez-Reche and E.Vives, Phys. Rev. B 67, 134421 (2003).
  * Ortín and Goicoechea (1998) J. Ortín and J. Goicoechea, Phys. Rev. B 58, 5628 (1998).
  * X.Illa et al. (2005) X.Illa, J.Ortín, and E.Vives, Phys. Rev. B 71, 184435 (2005).
  * Vilar and Rubi (2008) J. M. G. Vilar and J. M. Rubi, Phys. Rev. Lett. 100, 020601 (2008).
  * A.Imparato and L.Peliti (2007) A.Imparato and L.Peliti, Comptes Rendus Physique 8, 556 (2007).
  * J.Horowitz and C.Jarzynski (2008) J.Horowitz and C.Jarzynski, Phys. Rev. Lett. 101, 098901 (2008).
  * Peliti (2008) L. Peliti, Phys. Rev. Lett. 101, 098903 (2008).