Oppenheimer-Snyder collapse

This Jupyter/SageMath notebook is relative to the lectures Geometry and physics of black holes.

The computations make use of tools developed through the SageManifolds project.

In [1]:
version()
Out[1]:
'SageMath version 9.5.beta9, Release Date: 2021-12-23'
In [2]:
%display latex

Due to spherical symmetry, we a considering a 2-dimensional cut at $(\theta,\varphi) = \mathrm{const}$:

In [3]:
M = Manifold(2, 'M', structure='Lorentzian')
EF.<ti, r> = M.chart(r"ti:\tilde{t} r:[0,+oo)")
EF
Out[3]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\left(M,({\tilde{t}}, r)\right)\]
In [4]:
EF.coord_range()
Out[4]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}{\tilde{t}} :\ \left( -\infty, +\infty \right) ;\quad r :\ \left[ 0 , +\infty \right)\]

Interior

In [5]:
I = M.open_subset('I')

We choose $\chi_{\rm s} = \pi/4$:

In [6]:
S.<eta,chi> = I.chart(r"eta:[0,pi]:\eta chi:[0,pi/4]:\chi")
S
Out[6]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\left(I,({\eta}, {\chi})\right)\]
In [7]:
S.coord_range()
Out[7]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}{\eta} :\ \left[ 0 , \pi \right] ;\quad {\chi} :\ \left[ 0 , \frac{1}{4} \, \pi \right]\]
In [8]:
S.coord_bounds()
Out[8]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\left(\left(\left(0, \mathrm{True}\right), \left(\pi, \mathrm{True}\right)\right), \left(\left(0, \mathrm{True}\right), \left(\frac{1}{4} \, \pi, \mathrm{True}\right)\right)\right)\]
In [9]:
chis = S.coord_bounds()[1][1][0]
chis
Out[9]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{4} \, \pi\]

The initial areal radius of the star in units of $m$:

In [10]:
r0 = 2/sin(chis)^2 
r0
Out[10]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}4\]

Global evolution curves

$r_{\rm s}(\eta)$ in units of $r_0$:

In [11]:
rsurf(eta) = (1 + cos(eta))/2
rsurf
Out[11]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}{\eta} \ {\mapsto}\ \frac{1}{2} \, \cos\left({\eta}\right) + \frac{1}{2}\]

Matter proper time $\tau$ in units of $\sqrt{r_0^3/(8 m_0)}$:ยต

In [12]:
tau(eta) = (eta + sin(eta))/2
tau
Out[12]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}{\eta} \ {\mapsto}\ \frac{1}{2} \, {\eta} + \frac{1}{2} \, \sin\left({\eta}\right)\]

Proper energy density $\rho$ in units of $6 m / (\pi r_0^3)$:

In [13]:
rho(eta) = 1 / (1 + cos(eta))^3
rho
Out[13]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}{\eta} \ {\mapsto}\ \frac{1}{{\left(\cos\left({\eta}\right) + 1\right)}^{3}}\]
In [14]:
graph = plot(tau, (0, pi), color='purple', thickness=2, linestyle='--',
             legend_label=r'$\tau \;  \sqrt{2 m / r_0^3}$',
             axes_labels=[r'$\eta$', r'$r_{\rm s}$, $\tau$, $\rho$'],
             ticks=[pi/4, None], tick_formatter=[pi, None], fontsize=12,
             axes_labels_size=1.4,
             frame=True, axes=False, gridlines=True) \
        + plot(rsurf, (0, pi), color='red', thickness=2, 
               legend_label=r'$r_{\rm s}(\tau)/r_0 = a(\tau)\, \sin\chi_{\rm s} / r_0$') \
        + plot(rho, (0, 2), color='blue', thickness=2, 
               legend_label=r'$\rho(\tau) \, \pi r_0^3/(6 m)$')
graph.set_legend_options(handlelength=2, labelspacing=0.5,
                         font_size=11)
graph.set_aspect_ratio(0.75)
graph.show(ymax=3.2)
graph.save('lem_OS_rs_rho_eta.pdf', ymax=3.2)
In [15]:
graph = parametric_plot((tau(eta), eta), (eta, 0, pi), color='orange', 
                        linestyle='--', thickness=2, 
                        legend_label=r'$\eta$',
                        axes_labels=[r'$\tau \; \sqrt{2m/r_0^3}$', 
                                     r'$\eta$, $r_{\rm s}$, $\rho$'], 
                        frame=True, axes=False, gridlines=True) \
        + parametric_plot((tau(eta), rsurf(eta)), (eta, 0, pi), color='red',
                          thickness=2, 
                          legend_label=r'$r_{\rm s}(\tau)/r_0 = a(\tau)\, \sin\chi_{\rm s} / r_0$') \
        + parametric_plot((tau(eta), rho(eta)), (eta, 0, 2), color='blue', 
                          thickness=2, 
                          legend_label=r'$\rho(\tau) \, \pi r_0^3/(6 m)$')
graph.set_legend_options(handlelength=2, labelspacing=0.5,
                         font_size=11)
graph.set_aspect_ratio(0.4)
graph.show(ymax=3.2)
graph.save('lem_OS_rs_rho_tau.pdf', ymax=3.2)

$\eta$ as a function of $\tau/m$:

In [16]:
def etaf(tau):
    def ff(et):
        return (et + sin(et))/sin(chis)^3 - tau
    return find_root(ff, 0, pi)

Value of the matter proper time $\tau$ at the end of the collapse, in units of $m$:

In [17]:
tau_end = pi/sin(chis)^3
tau_end, n(tau_end)
Out[17]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\left(2 \, \sqrt{2} \pi, 8.88576587631673\right)\]

Check:

In [18]:
etaf(tau_end) == n(pi)
Out[18]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}\]
In [19]:
plot(etaf, 0, tau_end)
Out[19]:

Plot of $\rho$ at the black hole's birth as a function of $\chi_s$

In [20]:
graph = plot(3/(4*pi)*sin(x)^6/(1 - cos(3*x))^3, (x, 0, pi/2), 
             color='blue', thickness=2, 
             axes_labels=[r'$\chi_{\rm s}$', r'$\rho_{\rm hb} m^2$'],
             ticks=[pi/8, None], tick_formatter=[pi, None],
             frame=True, axes=False, gridlines=True)
graph.show(ymax=0.25)
graph.save("lem_OS_rho_hb.pdf", ymax=0.25)

Plot of the interior region in the $(\eta,\chi)$ plane

Constant $\chi$ lines

In [21]:
graph = S.plot(S, ambient_coords=(chi, eta), style={eta: '-', chi: '--'},
               number_values={eta: 2, chi: 9})

Constant $\tau$ lines

In [22]:
tau_sel = [numerical_approx(i*tau_end/8) for i in range(9)]
tau_sel
Out[22]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\left[0.000000000000000, 1.11072073453959, 2.22144146907918, 3.33216220361878, 4.44288293815837, 5.55360367269796, 6.66432440723755, 7.77504514177714, 8.88576587631673\right]\]
In [23]:
eta_sel = [etaf(tau0) for tau0 in tau_sel]
eta_sel
Out[23]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\left[0.0, 0.1969852777521357, 0.39790775618592655, 0.6073795752895496, 0.8317111935799412, 1.0810453707259005, 1.3752523669869274, 1.7683424316081915, 3.141592653589793\right]\]
In [24]:
iso_taus = [I.curve({S: [eta0, chi]}, (chi, 0, chis)) for eta0 in eta_sel]
In [25]:
for iso_tau in iso_taus:
    graph += iso_tau.plot(chart=S, ambient_coords=(chi, eta), color='red',
                          style='--')

The stellar surface:

In [26]:
surf = I.curve({S: [eta, chis]}, (eta, 0, pi)) 
graph += surf.plot(S, ambient_coords=(chi, eta), thickness=3)

The star at $\eta=0$:

In [27]:
init= I.curve({S: [0, chi]}, (chi, 0, chis)) 
graph += init.plot(S, ambient_coords=(chi, eta), 
                   color='magenta', thickness=3)
graph.show(figsize=8)

Null radial geodesics in the interior

Functions initializing the internal null geodesics

In [28]:
def geod_int_out(etas):
    r"""
    Internal outgoing radial null geodesic ending at
    (eta, chi) = (etas, chis)
    """
    chi_min = 0 if etas > chis else chis - etas
    return I.curve({S: [chi - chis + etas, chi]}, (chi, chi_min, chis))

def geod_int_in(etas):
    r"""
    Internal ingoing radial null geodesic starting at
    (eta, chi) = (etas, chis)
    """
    chi_min = 0 if etas < pi - chis else etas - pi + chis
    return I.curve({S: [chis - chi + etas, chi]}, (chi, chi_min, chis))

The event horizon

In [29]:
hor_int = geod_int_out(pi - 2*chis)
graph += hor_int.plot(chart=S, ambient_coords=(chi, eta),
                      color='black', thickness=3)

The final singularity

In [30]:
matter_sing = I.curve({S: [pi, chi]}, (chi, 0, chis))

graph += line([(0 + i*chis/8, pi) for i in range(9)], 
              thickness=3, color='darkorange',
              marker='D',  markersize=10)
graph.show(frame=True, axes_pad=0, figsize=8)

Eddington-Finkelstein-type coordinates for the interior

In [31]:
EFI = EF.restrict(I)
I.atlas()
Out[31]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left(I,({\eta}, {\chi})\right), \left(I,({\tilde{t}}, r)\right)\right]\]
In [32]:
S_to_EF = S.transition_map(EFI, 
                           [2*(sqrt(r0/2 - 1)*(eta + r0/4*(eta + sin(eta))) 
                             + 2*ln(cos(eta/2) + sin(eta/2)/sqrt(r0/2 - 1))),
                            r0*sin(chi)/sin(chis)/2*(1 + cos(eta))])
S_to_EF.display()
Out[32]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\left\{\begin{array}{lcl} {\tilde{t}} & = & 4 \, {\eta} + 4 \, \log\left(\cos\left(\frac{1}{2} \, {\eta}\right) + \sin\left(\frac{1}{2} \, {\eta}\right)\right) + 2 \, \sin\left({\eta}\right) \\ r & = & 2 \, \sqrt{2} {\left(\cos\left({\eta}\right) + 1\right)} \sin\left({\chi}\right) \end{array}\right.\]
In [33]:
S_to_EF(eta, chi)
Out[33]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\left(4 \, {\eta} + 4 \, \log\left(\cos\left(\frac{1}{2} \, {\eta}\right) + \sin\left(\frac{1}{2} \, {\eta}\right)\right) + 2 \, \sin\left({\eta}\right), 2 \, \sqrt{2} \cos\left({\eta}\right) \sin\left({\chi}\right) + 2 \, \sqrt{2} \sin\left({\chi}\right)\right)\]
In [34]:
rr(eta, chi) = S_to_EF(eta, chi)[1]
rr
Out[34]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\left( {\eta}, {\chi} \right) \ {\mapsto} \ 2 \, \sqrt{2} \cos\left({\eta}\right) \sin\left({\chi}\right) + 2 \, \sqrt{2} \sin\left({\chi}\right)\]

Adding isocontours $r = \mathrm{const}$:

In [35]:
graph += contour_plot(rr(eta, chi), (chi, 0, chis), (eta, 0, 3.13), 
                      contours=[0.5, 1., 1.5, 2., 2.5, 3., 3.5], 
                      fill=False, labels=True, label_colors='black',
                      label_inline=True, label_inline_spacing=10)
In [36]:
graph.show(frame=True, axes_pad=0, ticks=[pi/8,pi/8], tick_formatter=(pi,pi), 
           fontsize=16, axes_labels_size=1.2, figsize=10)
In [37]:
ITH = I.curve({S: [pi - 2*chi, chi]}, (chi, 0, chis))

graph += ITH.plot(chart=S, ambient_coords=(chi, eta), color='blue', 
                  thickness=3, style=':', label_axes=False)
In [38]:
graph_trap = copy(graph)

A selection of null geodesics:

In [39]:
geod_int_out_sel = [geod_int_out(etas) for etas in eta_sel]
geod_int_in_sel = [geod_int_in(etas) for etas in eta_sel]

for geod in geod_int_out_sel:
    graph += geod.plot(chart=S, ambient_coords=(chi, eta),
                       color='green', thickness=2)

for geod in geod_int_in_sel[:-1]:
    graph += geod.plot(chart=S, ambient_coords=(chi, eta),
                       color='green',  thickness=1, style='--')

graph.show(frame=True, axes=False, axes_pad=0, ticks=[pi/8,pi/4], 
           tick_formatter=(pi,pi), fontsize=16, axes_labels_size=1.2, 
           xmax=pi/4+0.015, ymin=-0.015, figsize=10)
graph.save('lem_OS_diag_int.pdf', frame=True, axes=False, axes_pad=0, 
           ticks=[pi/8,pi/4], tick_formatter=(pi,pi), fontsize=16, 
           axes_labels_size=1.2, xmax=pi/4+0.015, ymin=-0.015, figsize=10)

Trapped surfaces

In [40]:
for i in range(9):
    etas = pi/8*i 
    graph_trap += geod_int_out(etas).plot(chart=S, 
                                          ambient_coords=(chi, eta),
                                          color='green', thickness=2)
    graph_trap += geod_int_in(etas).plot(chart=S, 
                                         ambient_coords=(chi, eta),
                                         color='green',  thickness=1, 
                                         style='--')
graph_trap.show(frame=True, axes_pad=0, ticks=[pi/8,pi/4], tick_formatter=(pi,pi), 
                fontsize=16, axes_labels_size=1.2, figsize=10)

Adding $r$ isocontours near the trapped surface at $(\eta,\chi) = \left(\frac{11\pi}{6}, \frac{3\pi}{16} \right)$:

In [41]:
tsp = I((11*pi/16, 3*pi/16), chart=S)
In [42]:
S(tsp)
Out[42]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\left(\frac{11}{16} \, \pi, \frac{3}{16} \, \pi\right)\]
In [43]:
EF(tsp)
Out[43]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\left(\frac{11}{4} \, \pi + 4 \, \log\left(\cos\left(\frac{11}{32} \, \pi\right) + \sin\left(\frac{11}{32} \, \pi\right)\right) + 2 \, \sin\left(\frac{5}{16} \, \pi\right), -2 \, \sqrt{2} \cos\left(\frac{5}{16} \, \pi\right) \sin\left(\frac{3}{16} \, \pi\right) + 2 \, \sqrt{2} \sin\left(\frac{3}{16} \, \pi\right)\right)\]
In [44]:
rts = n(EF(tsp)[1])
rts
Out[44]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}0.698372454547306\]
In [45]:
graph_trap += contour_plot(rr(eta, chi), (chi, 0, chis), (eta, 0, 3.13), 
                           contours=[rts - 0.05, rts, rts + 0.05], 
                           fill=False, labels=True, label_colors='black',
                           label_inline=True, label_inline_spacing=10)
graph_trap += tsp.plot(chart=S, ambient_coords=(chi, eta), color='brown',
                       label=' ', size=60) \
              + text(r"$\mathscr{S}$", (0.63, 2.16), color='brown', fontsize=18)
graph_trap.show(frame=True, axes_pad=0, ticks=[pi/8,pi/8], tick_formatter=(pi,pi),
                figsize=8, xmax=pi/4+0.005, ymin=1.4, ymax=2.5)
graph_trap.save('lem_OS_diag_int_zoom.pdf', frame=True, axes_pad=0, ticks=[pi/8,pi/8], 
                tick_formatter=(pi,pi), figsize=8, xmax=pi/4+0.005, ymin=1.4, ymax=2.5)

Plots in terms of Eddington-Finkelstein coordinates

In [46]:
graph = S.plot(EFI, ambient_coords=(r, ti), style={eta: '-', chi: '--'},
               number_values={eta: 2, chi: 9}, label_axes=False) 

# The stellar surface:
graph += surf.plot(EFI, ambient_coords=(r, ti), fixed_coords={chi: chis}, 
                   thickness=3, label_axes=False)

# Constant tau lines:
for iso_tau in iso_taus:
    graph += iso_tau.plot(chart=EFI, ambient_coords=(r, ti), color='red',
                          style='--', label_axes=False)

# The initial star:
graph += init.plot(chart=EFI, ambient_coords=(r, ti), 
                   fixed_coords={chi: chis}, color='magenta', 
                   thickness=3, label_axes=False)
    
graph.axes_labels([r'$r/m$', r'$\tilde{t}/m$'])
graph.show(figsize=10)

Coordinate $\tilde{t}$ of the collapse end point

The collapse end point correspond to $\eta = \pi$; it is thus given by I((pi, 0)) and its $\tilde{t}$ coordinate is obtained via the chart EF:

In [47]:
EF(I((pi, 0)))
Out[47]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\left(4 \, \pi, 0\right)\]
In [48]:
ti_end = EF(I((pi, 0)))[0]
ti_end
Out[48]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}4 \, \pi\]

Drawing the singularity

Maximal value of $\tilde{t}$ for plots:

In [49]:
timax = 17
In [50]:
graph += line([(0, ti_end + i*(40 - ti_end)/80) for i in range(81)], 
              thickness=3, color='darkorange', marker='D', markersize=10)
graph.set_axes_range(xmin=0, xmax=10, ymin=0, ymax=timax)
graph.axes_labels([r'$r/m$', r'$\tilde{t}/m$'])
graph.show(frame=True, figsize=10, axes_pad=0)

The event horizon

In [51]:
graph += hor_int.plot(chart=EFI, ambient_coords=(r, ti), color='black', 
                      thickness=3, label_axes=False)
In [52]:
hor_int(chis)
Out[52]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\mbox{Point on the 2-dimensional Lorentzian manifold M}\]
In [53]:
EF(hor_int(chis))
Out[53]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\left(2 \, \pi + 2 \, \log\left(2\right) + 2, 2\right)\]
In [54]:
tiH = EF(hor_int(chis))[0]
tiH, n(tiH)
Out[54]:
\[\newcommand{\Bold}[1]{\mathbf{#1}}\left(2 \, \pi + 2 \, \log\left(2\right) + 2, 9.66947966829948\right)\]
In [55]:
hor_ext = M.curve({EF: (ti, 2)}, (ti, tiH, 40))
In [56]:
graph += hor_ext.plot(chart=EF, ambient_coords=(r, ti), color='black', 
                      thickness=3, label_axes=False)

Adding the inner trapping horizon:

In [57]:
graph += ITH.plot(EFI, ambient_coords=(r, ti), color='blue', 
                  thickness=3, style=':', label_axes=False)
graph.set_axes_range(xmin=0, xmax=10, ymin=0, ymax=timax)
graph.show(figsize=10)
In [58]:
graph_delay = copy(graph)

Functions initializing the external null geodesics

In [59]:
def geod_ext_out_I(tis, rs, rm):
    r"""
    External outgoing radial null geodesic starting from
    (ti, r) = (tis, rs) in region I
    """
    return M.curve({EF: [tis + r - rs + 4*ln(abs((r - 2)/(rs - 2))), r]}, (r, rs, rm))

lamb = var('lamb', latex_name=r'\lambda')
def geod_ext_out_II(tis, rs, rm):
    r"""
    External outgoing radial null geodesic starting from
    (ti, r) = (tis, rs) in region II
    """
    return M.curve({EF: [tis - lamb - rs + 4*ln(abs((-lamb - 2)/(rs - 2))), -lamb]}, (lamb, -rs, -rm))

def geod_ext_in(tis, rs, rm):
    r"""
    External ingoing radial null geodesic arriving at
    (ti, r) = (tis, rs)
    """
    return M.curve({EF: [tis - r + rs, r]}, (r, rs, rm))
In [60]:
def plot_geod_out(etas, chart=EF, plot_tangent=True, plot_int=True,
                  color='green', thickness=1.5):
    r"""
    Plot of an outgoing radial null geodesic
    """
    # Internal part of the geodesic
    gint = geod_int_out(etas)
    # External part
    tis, rs = EF(gint(chis))
    if rs < 2:
        rm = 0
        gext = geod_ext_out_II(tis, rs, rm)
        pv = -0.5*rs
    else:
        rm = 20 
        gext = geod_ext_out_I(tis, rs, rm)
        pv = 1.4*rs
    ambc = (chart[1], chart[0])
    graph = gext.plot(chart, ambient_coords=ambc, color=color, 
                      thickness=thickness, label_axes=False)
    if plot_int:
        graph += gint.plot(chart.restrict(I), ambient_coords=ambc, color=color, 
                           thickness=thickness, label_axes=False) 
    if plot_tangent:
        try:
            vtan = gext.tangent_vector_field().at(gext.domain()(n(pv)))
            graph += vtan.plot(chart, ambient_coords=ambc, color=color, 
                               scale=0.01, arrowsize=4)
        except ValueError:
            pass
    return graph
In [61]:
for etas in eta_sel:
    graph += plot_geod_out(etas)
In [62]:
graph.axes_labels([r'$r/m$', r'$\tilde{t}/m$'])
graph.set_axes_range(xmin=0, xmax=10, ymin=0, ymax=timax)
graph.show(figsize=10)
In [63]:
def plot_geod_in(etas, chart=EF, color='green', thickness=1):
    r"""
    Plot of an ingoing radial null geodesic
    """
    # Internal part of the geodesic
    gint = geod_int_in(etas)
    # External part
    tis, rs = EF(gint(chis))
    rm = 20
    gext = geod_ext_in(tis, rs, rm)
    ambc = (chart[1], chart[0])
    return gint.plot(chart.restrict(I), ambient_coords=ambc, 
                     color=color, thickness=thickness, style='--', 
                     label_axes=False) \
           + gext.plot(chart, ambient_coords=ambc, color=color, 
                       thickness=thickness, style='--', 
                       label_axes=False)
In [64]:
for etas in eta_sel:
    graph += plot_geod_in(etas)

for ti0 in [14, 18, 22, 26]:
    graph += geod_ext_in(ti0, 0, 20).plot(EF, ambient_coords=(r, ti), 
                                          color='green', thickness=1, 
                                          style='--', label_axes=False)
graph_zoom = copy(graph)
graph += text(r'$\mathscr{H}$', (2.5, 16.4), fontsize=20, color='black')
graph.axes_labels([r'$r/m$', r'$\tilde{t}/m$'])
graph.set_axes_range(xmin=0, xmax=10, ymin=-0.03, ymax=timax)
graph.show(frame=True, axes=False, axes_pad=0, figsize=10)
graph.save('lem_OS_diag_EF.pdf', frame=True, axes=False, 
           axes_pad=0, figsize=10)
In [65]:
graph_zoom += plot_geod_out(3*pi/4) + plot_geod_in(5*pi/8)
graph_zoom += tsp.plot(chart=EFI, ambient_coords=(r, ti), color='brown',
                       label=' ', size=50) \
              + text(r"$\mathscr{S}$", (0.8, 11.6), color='brown', fontsize=18)
graph_zoom += line([(0, ti_end + i*(15 - ti_end)/20) for i in range(21)], 
                   thickness=3, color='darkorange', marker='D', markersize=10)
In [66]:
graph_zoom.axes_labels([r'$r/m$', r'$\tilde{t}/m$'])
graph_zoom.show(frame=True, figsize=8, axes_pad=0,
                xmin=0, xmax=3, ymin=9, ymax=13)
graph_zoom.save('lem_OS_diag_EF_zoom.pdf',frame=True, 
                figsize=8, axes_pad=0, xmin=0, xmax=3, 
                ymin=9, ymax=13)