Sampling flux vacua

What’s in this notebook? This notebook is the practical companion to NB06 — ISD sampling principle. It shows how to turn ISD-completed flux seeds into refined SUSY vacua, how to use the high-level sampling wrappers, and how to run a compact non-SUSY critical-point search.

In this notebook, you will learn:

• How to refine integer-rounded ISD fluxes to true SUSY vacua with linearised_shifts and Newton refinement.

• How to run the same SUSY workflow through sample_SUSY_flux_vacua.

• How the non-SUSY sample_critical_points workflow differs from the SUSY F-term solve.

• How to inspect sampled ensembles with lightweight comparison plots.

Non-goals: detailed bounded-flux enumeration, solver benchmarking, and production performance tuning live in NB08 — Flux bounding and NB16 — Sampling benchmarks. Literature-data comparisons live in NB17 — Flux bounding case studies.

Prerequisites: NB04 — Sampling module, NB05 — Finding flux vacua, and NB06 — ISD sampling principle.

Recommended paths through this notebook

This tutorial is now scoped to the sampling workflows themselves. If you only need a working SUSY sampler, read Workflows A and B. If you want to see how the same ISD-biased idea feeds a scalar-potential critical-point search, add Workflow C.

If you want to…	Read
inspect the ISD-to-Newton pipeline	Workflow A
use the shortest production-style SUSY API	Workflow B
sample non-SUSY critical points	Workflow C
compare solvers or tuning choices	NB16
understand bounded enumeration	NB08

Outline

1. Setup

2. Workflow A: SUSY ISD sampling from input fluxes

3. Workflow B: SUSY ISD sampling via the wrapper API

4. Workflow C: Non-SUSY ISD sampling

5. Take-aways

6. Further reading

Setup

# General imports
import warnings, sys, time
import numpy as np
from tqdm.auto import tqdm
from scipy.optimize import root

# JAX imports
import jax
import jax.numpy as jnp
jax.config.update("jax_enable_x64", True)

# Plotting
import seaborn as sn
import matplotlib.pyplot as plt
cmap = sn.color_palette("viridis", as_cmap=True)

# JAXVacua
import jaxvacua as jvc
from stringjax_tools import vmapping_func

sys.path.append("../")
from plot_utils import make_overview_plots

warnings.filterwarnings("ignore")
# Keep the default notebook execution short. The h12=3 examples are useful
# scaling checks, but compile/run longer and are disabled by default.
RUN_H12_3_EXAMPLES = False

Canonical fixture. Throughout this notebook we use the degree-18 hypersurface in \(\mathbb{CP}^{1,1,1,6,9}\) (\(h^{2,1}=2\), \(Q_{D3}=276\), model_ID=1) as the running example. The same fixture is reused in every workflow and benchmark so that results are directly comparable. The h12=3 mini-examples at the end of §A and §C illustrate that the workflows generalise without code changes.

h12 = 2
model = jvc.FluxVacuaFinder(
    h12=h12,
    Q=276,
    model_ID=1,
    maximum_degree=2,
    map_to_fd=True,   # map returned tau to the SL(2,Z) fundamental domain
)

We build a data_sampler with mildly restrictive flux and moduli bounds that we reuse across the practical workflows below. The construction is identical to NB04: the sampler controls where initial moduli, axio-dilaton values and flux seeds are drawn; it does not by itself solve the F-term equations.

sampler = jvc.data_sampler(
    model,
    flux_bounds=[-5, 5],
    moduli_bounds=[0, 10],
    axion_bounds=[-0.5, 0.5],
)

For convenient ensemble evaluation we vectorise the F-term function model.DW once and reuse it throughout. The cross-link to the ISD principle: the algorithmic intro to ISD sampling lives in NB06; here we focus on the workflow of producing vacua at scale.

objective_fct = model.DW

Workflow A: SUSY ISD sampling from input fluxes

Workflow A starts from a fixed integer flux choice and searches for all moduli configurations that yield a SUSY vacuum compatible with it. This is the complement of Workflow B (moduli-first); together they cover the two natural ways of populating an ensemble. The algorithmic primitive is the same in both — the ISD-rounded shift from NB06 ISD principle — but the order in which fluxes and moduli are sampled differs.

The algorithm

The pipeline has two stages, each vmap’d twice (once over moduli starting points and once over candidate flux vectors):

1. find_solution_init_vmap — apply model.linearised_shifts to round the input flux into an ISD-compatible integer flux and produce near-vacuum starting moduli.

2. find_solution_steps_vmap — iterate model.linearised_shifts to drive the residual \(|D_I W|\) to zero (Newton-like).

Below we wire up both stages, run a few iterations on two flux samples, and confirm the residual shrinks at each step.

# Common kwargs for the optimiser. "Fflux" fixes the RR-fluxes as input;
# the linearised_shifts call then solves for the NS-NS half.
mode = "Fflux"
kwargs = {"mode": mode, "Q": model.Q(),
          "constraints": None, "remove_NANs": True}

# Double vmap of model.DW — once over moduli starting points, once over fluxes.
DW_vv = jax.vmap(jax.vmap(model.DW))

The two doubly-vmapped optimisers — one for the ISD initialisation step (find_solution_init_vmap) and one for the Newton refinement step (find_solution_steps_vmap) — are built with the vmapping_func helper imported in §1:

# Step 1: ISD initialisation. Outer vmap: moduli starting points. Inner: same flux.
find_solution_init = vmapping_func(
    model.linearised_shifts, in_axes=(0, 0, None), return_flag=False, **kwargs)
find_solution_init_vmap = vmapping_func(
    find_solution_init, in_axes=(None, None, 0))

# Step 2: Newton refinement. Both axes vmapped.
find_solution_steps = vmapping_func(
    model.linearised_shifts, in_axes=(0, 0, 0), return_flag=True, **kwargs)
find_solution_steps_vmap = vmapping_func(
    find_solution_steps, in_axes=(0, 0, 0))

We supply two starting points and two candidate flux vectors. The doubly-vmapped call evaluates every (moduli, flux) pair in parallel.

moduli0 = np.array([
    [ 0.22924045 + 3.00550378j, -0.46016314 + 2.28651085j],
    [-0.04468382 + 0.63834065j,  0.25839392 + 2.36493542j],
])
tau0    = np.array([0.18254317 + 6.81361114j, 0.49828078 + 8.33716844j])
fluxes0 = np.array([[ 1,  5,  1,  2,  4, -1],
                    [ 5, -3,  0, -4, -4, -3]])

# Step 1: ISD initialisation rounds the fluxes and shifts the moduli.
moduli, tau, fluxes = find_solution_init_vmap(moduli0, tau0, fluxes0)

Iterating the Newton-like step find_solution_steps_vmap then drives \(|D_I W|\) toward zero. We track the per-flux residual at each step:

moduli1, tau1, fluxes1, _ = find_solution_steps_vmap(moduli, tau, fluxes)

for i in range(10):
    moduli1, tau1, fluxes1, _ = find_solution_steps_vmap(moduli1, tau1, fluxes1)
    DW = DW_vv(moduli1, jnp.conj(moduli1), tau1, jnp.conj(tau1), fluxes1)
    res = np.sum(np.abs(DW), axis=2)
    print(f"step {i:2d}   |DW| flux1: {res[0]}   flux2: {res[1]}")

Batch workflow

The same two-stage ISD initialisation + Newton refinement pipeline can be batched over many starting points. The default cell below uses a deliberately small target so that the notebook remains executable as documentation. Increase N, vmap_dim_pts, and the number of initial fluxes only when you are preparing a production run.

rns_key = jax.random.PRNGKey(42)

moduli, tau, fluxes, residual = model.sample_SUSY_vacua_from_fluxes(
    fluxes_init=fluxes0,
    N=50,
    sampler=sampler,
    rns_key=rns_key,
    moduli_sampling_mode="cone",
    max_iters=10,
    objective_fct=DW_vv,
    optimiser_init=find_solution_init_vmap,
    optimiser_steps=find_solution_steps_vmap,
    vmap_dim_pts=25,
    mode=mode,
    print_progress=True,
)

DW = model.DW(moduli, jnp.conj(moduli), tau, jnp.conj(tau), fluxes)
print(f"all |DW| residuals below 1e-10: "
      f"{np.all(np.sum(np.abs(DW), axis=1) < 1e-10)}")

The fixed-flux run returns a small ensemble of refined vacua. A quick overview plot is often the fastest way to catch obvious pathologies: moduli outside the intended patch, strange axio-dilaton values, or a tadpole distribution dominated by a few outliers.

make_overview_plots(
    model, moduli, tau, fluxes,
    moduli_range=[[-0.5, 0.5], [1, 4]],
    tau_range=[[-0.5, 0.5], [3, 10]],
    W0_range=[-10, 10],
    use_normal_w0=False,
)

Improved constrained sampling — input fluxes across a moduli region

By passing fluxes_init=None we instead let the sampler draw fluxes internally, while still constraining the moduli to a chosen sub-region of the Kähler cone. This is the flux-first analog of moduli-first sampling: instead of fixing fluxes and searching moduli, we fix a region of moduli space and let the sampler explore fluxes consistent with it.

moduli, tau, fluxes, residuals = model.sample_SUSY_vacua_from_fluxes(
    fluxes_init=None,
    N=30,
    sampler=sampler,
    rns_key=rns_key,
    max_iters=10,
    moduli_sampling_mode="cone",
    Q=276,
    objective_fct=DW_vv,
    optimiser_init=find_solution_init_vmap,
    optimiser_steps=find_solution_steps_vmap,
    vmap_dim_flux=50,
    vmap_dim_pts=25,
    mode=mode,
    print_progress=True,
)

The constrained run uses internally sampled fluxes rather than fixed input fluxes. The same diagnostic plot now checks the output of the region-restricted search.

make_overview_plots(
    model, moduli, tau, fluxes,
    moduli_range=[[-0.5, 0.5], [1, 4]],
    tau_range=[[-0.5, 0.5], [3, 10]],
    W0_range=[-10, 10],
    use_normal_w0=False,
)

Higher-\(h^{1,2}\) example: SUSY workflow at \(h^{1,2}=3\)

The code path is the same in higher-dimensional moduli spaces, but compilation and sampling costs rise quickly. This example is therefore guarded by RUN_H12_3_EXAMPLES; leave it off during routine notebook execution.

if RUN_H12_3_EXAMPLES:
    # h12=3 example (model_ID=4 has hyperplanes + sensible sampling region).
    model_3 = jvc.FluxVacuaFinder(h12=3, model_ID=4, maximum_degree=2)

    sampler_3 = jvc.data_sampler(
        model_3,
        flux_bounds=[-5, 5],
        moduli_bounds=(2., 8.),
        dilaton_bounds=(np.sqrt(3) / 2, 10.),
        axion_bounds=(-0.5, 0.5),
        seed=42,
    )

    # Rebuild the optimisers with the new model.
    kwargs_3 = {**kwargs, "Q": model_3.Q()}
    _init_3 = vmapping_func(model_3.linearised_shifts, in_axes=(0, 0, None),
                           return_flag=False, **kwargs_3)
    init_3_vmap = vmapping_func(_init_3, in_axes=(None, None, 0))
    _step_3 = vmapping_func(model_3.linearised_shifts, in_axes=(0, 0, 0),
                           return_flag=True, **kwargs_3)
    step_3_vmap = vmapping_func(_step_3, in_axes=(0, 0, 0))
    DW_vv_3 = jax.vmap(jax.vmap(model_3.DW))

    moduli3, tau3, fluxes3, res3 = model_3.sample_SUSY_vacua_from_fluxes(
        fluxes_init=None,
        N=50,
        sampler=sampler_3,
        rns_key=jax.random.PRNGKey(42),
        max_iters=10,
        moduli_sampling_mode="cone",
        Q=276,
        objective_fct=DW_vv_3,
        optimiser_init=init_3_vmap,
        optimiser_steps=step_3_vmap,
        vmap_dim_flux=25,
        vmap_dim_pts=10,
        mode=mode,
        print_progress=True,
    )
    DW3 = model_3.DW(moduli3, jnp.conj(moduli3),
                              tau3, jnp.conj(tau3), fluxes3)
    print(f"h12=3 — {len(moduli3)} vacua; "
          f"all |DW| < 1e-10: {np.all(np.sum(np.abs(DW3), axis=1) < 1e-10)}")
else:
    print("RUN_H12_3_EXAMPLES = False — skipping the h12=3 SUSY example.")

Workflow B: SUSY ISD sampling via the wrapper API

What’s new vs §A. §A wired the two-stage pipeline manually: build two doubly-vmapped optimisers (find_solution_init_vmap and find_solution_steps_vmap), call them in sequence, iterate. §B wraps the same pipeline behind a single high-level driver, model.sample_SUSY_flux_vacua, that handles the moduli draws, ISD rounding, and Newton refinement in one call. The trade-off: less fine-grained control over the inner loop, but a much shorter user-facing snippet that suffices for routine production runs.

ISD sampling with a single optimiser

sample_SUSY_flux_vacua wraps the manual pipeline from Workflow A behind one API call. The default target is small enough for a tutorial run; use larger N and vmap_dim values for production sampling.

rns_key = jax.random.PRNGKey(42)

kwargs = {"Q": model.Q(), "return_flag": True,
          "constraints": None, "remove_NANs": True,
          "in_axes": (0, 0, 0)}

linearised_shifts_F = vmapping_func(
    model.linearised_shifts, mode="Fflux", **kwargs)

moduli, tau, fluxes, residuals = model.sample_SUSY_flux_vacua(
    N=50,
    rns_key=rns_key,
    sampler=sampler,
    optimiser=linearised_shifts_F,
    objective_fct=objective_fct,
    vmap_dim=50,
    print_progress=True,
)

The wrapper produces the same kind of output as Workflow A, so we use the same compact overview plot to inspect the returned ensemble.

make_overview_plots(model, moduli, tau, fluxes,
                        moduli_range=[[-0.5, 0.5], [1, 4]],
                        tau_range=[[-0.5, 0.5], [3, 10]],
                        W0_range=[-10, 10],
                        use_normal_w0=False)

Running with multiple optimisers

Different ISD-completion modes can be tried in sequence. This is helpful when one mode has low yield in the chosen region, but each extra optimiser adds compilation and evaluation cost.

linearised_shifts_ISD = vmapping_func(
    model.linearised_shifts, mode="ISD",   **kwargs)
linearised_shifts_H   = vmapping_func(
    model.linearised_shifts, mode="Hflux", **kwargs)
linearised_shifts_F   = vmapping_func(
    model.linearised_shifts, mode="Fflux", **kwargs)

optimisers = [linearised_shifts_ISD, linearised_shifts_H, linearised_shifts_F]

moduli, tau, fluxes, residuals = model.sample_SUSY_flux_vacua(
    N=30,
    rns_key=rns_key,
    sampler=sampler,
    optimisers=optimisers,
    objective_fct=objective_fct,
    vmap_dim=50,
    print_progress=True,
)

Pooling several optimisers changes the sampled flux distribution. Plotting the output makes the comparison concrete before doing any statistical analysis.

make_overview_plots(model, moduli, tau, fluxes,
                        moduli_range=[[-0.5, 0.5], [1, 4]],
                        tau_range=[[-0.5, 0.5], [3, 10]],
                        W0_range=[-10, 10],
                        use_normal_w0=False)

Constrained sampling — moduli-region restrictions

Constraints let the sampler reject solutions outside a target moduli region. The example keeps the candidate count small and is meant to demonstrate the API, not to estimate final yields.

@jit
def constraints_model(moduli, tau, flux):
    return (jnp.all(jnp.imag(moduli) < 3., axis=0)
            & jnp.all(jnp.imag(moduli) > 2., axis=0)
            & (jnp.imag(tau) > 4))

kwargs_C = {"Q": model.Q(), "return_flag": True,
            "constraints": constraints_model, "remove_NANs": True,
            "in_axes": (0, 0, 0)}
linearised_shifts_F_v = vmapping_func(
    model.linearised_shifts, mode="Fflux", **kwargs_C)

# Tighter sampler matching the constraint region.
sampler_constrained = jvc.data_sampler(
    model, flux_bounds=[-5, 5],
    moduli_bounds=[2, 3], dilaton_bounds=[4, 10])

moduli, tau, fluxes, residuals = model.sample_SUSY_flux_vacua(
    N=30,
    rns_key=rns_key,
    sampler=sampler_constrained,
    optimiser=linearised_shifts_F_v,
    objective_fct=objective_fct,
    vmap_dim=50,
    print_progress=True,
)

The constrained sampler should visibly concentrate the ensemble in the requested region. The plot below is a quick sanity check of the imposed bounds.

make_overview_plots(model, moduli, tau, fluxes,
                        moduli_range=[[-0.5, 0.5], [1, 4]],
                        tau_range=[[-0.5, 0.5], [3, 10]],
                        W0_range=[-10, 10],
                        use_normal_w0=False)

Workflow C: Non-SUSY ISD sampling

What’s new vs §A/§B. Workflows A and B solve the F-term conditions \(D_I W = 0\) (SUSY vacua) using FluxVacuaFinder’s Newton + ISD pipeline. Workflow C targets the broader class of critical points of the scalar potential, \(\partial_I V = 0\), which includes non-SUSY vacua where \(D_I W \neq 0\). The non-SUSY workflow lives directly on FluxVacuaFinder via :func:sample_critical_points; the inner solve targets \(\partial V = 0\) instead of \(DW = 0\).

Background — non-SUSY critical points

SUSY flux vacua satisfy \(D_I W = 0\) (F-flatness), which implies \(\partial_I V = 0\) automatically. But \(V\) also has critical points where \(\partial_I V = 0\) while \(D_I W \neq 0\) — these are non-SUSY vacua. Two scalar-potential variants are commonly used:

• No-scale \(V_{\rm ns} = e^K |DW|^2 \geq 0\): all critical points are either SUSY minima (\(DW = 0\), \(V = 0\)) or non-SUSY critical points (\(DW \neq 0\), \(V > 0\)). Empirically ~85% of the critical points found are minima (not saddles).

• Full SUGRA \(V = e^K(|DW|^2 - 3|W|^2)\): can be negative, producing AdS minima and more saddles.

sample_critical_points defaults to the no-scale potential (noscale=True); switch to the full SUGRA potential via noscale=False.

Quick start: finding non-SUSY critical points

This is the scalar-potential analogue of the SUSY examples above: the sampler proposes ISD-biased flux candidates, but the solver targets \(\partial_I V = 0\) rather than \(D_I W = 0\). The default target is intentionally small; use NB16 for solver and throughput benchmarks.

sampler_nonSUSY = jvc.data_sampler(
    model,
    moduli_bounds=(2., 5.),
    dilaton_bounds=(np.sqrt(3) / 2, 10.),
    axion_bounds=(-0.5, 0.5),
    seed=42,
)

Nmax = 200  # tadpole bound for sampled fluxes

# moduli_max=100 (default) filters runaway solutions with |z_i| > 100.
finder = model

results = finder.sample_critical_points(
    Q=Nmax, noscale=True, sampler=sampler_nonSUSY, n_target=10,
    n_batch=50,
    max_batches=3,
    isd_mode="ISD-",
    solver="scipy",
    verbose=True,
)

Analysing the results

Each result is a dict with keys moduli, tau, flux, V, |DW|, is_susy, is_minimum, Nflux, etc. The lines below classify the batch into SUSY-vs-non-SUSY and minima-vs-saddles, then print the first ten entries.

n_susy = sum(1 for r in results if r['is_susy'])
n_nonsusy = len(results) - n_susy
n_min = sum(1 for r in results if r['is_minimum'])
n_sad = len(results) - n_min

print(f"Total critical points: {len(results)}")
print(f"  SUSY:     {n_susy} "
      f"({100*n_susy/max(len(results),1):.0f}%)")
print(f"  non-SUSY: {n_nonsusy} "
      f"({100*n_nonsusy/max(len(results),1):.0f}%)")
print(f"  minima:   {n_min} "
      f"({100*n_min/max(len(results),1):.0f}%)")
print(f"  saddle:   {n_sad} "
      f"({100*n_sad/max(len(results),1):.0f}%)")

if results:
    print(f"\n{'#':>3} {'V':>12} {'|DW|':>12} {'min_eig':>10} "
          f"{'type':>12} {'Nflux':>6}")
    print("-" * 60)
    for i, r in enumerate(results[:10]):
        label = "SUSY" if r['is_susy'] else "non-SUSY"
        kind = "minimum" if r['is_minimum'] else "saddle"
        print(f"{i:>3} {r['V']:>12.3e} {r['|DW|']:>12.3e} "
              f"{r.get('min_eig', float('nan')):>10.2e} "
              f"{label+'/'+kind:>12} "
              f"{r.get('Nflux', 0):>6.0f}")

For non-SUSY searches the most useful first diagnostic is not a moduli-space scatter, but the distribution of residuals, vacuum energies and tadpoles. The following small plot is intentionally cheap and works directly from the result dictionaries returned by sample_critical_points.

if results:
    V_vals = np.array([float(r["V"]) for r in results])
    DW_vals = np.array([float(r["|DW|"]) for r in results])
    N_vals = np.array([float(r.get("Nflux", np.nan)) for r in results])

    fig, axes = plt.subplots(1, 3, figsize=(11, 3.2), dpi=130)
    axes[0].hist(np.log10(np.maximum(DW_vals, 1e-16)), bins=20, color="steelblue", edgecolor="white")
    axes[0].set_xlabel(r"$\log_{10}|DW|$")
    axes[0].set_ylabel("count")
    axes[0].set_title("F-term residual")

    axes[1].hist(V_vals, bins=20, color="coral", edgecolor="white")
    axes[1].set_xlabel(r"$V$")
    axes[1].set_title("Vacuum energy")

    axes[2].hist(N_vals[np.isfinite(N_vals)], bins=20, color="mediumseagreen", edgecolor="white")
    axes[2].set_xlabel(r"$N_{\rm flux}$")
    axes[2].set_title("Tadpole")

    plt.tight_layout()
    plt.show()
else:
    print("No critical points returned; increase n_batch or max_batches for a denser run.")

Higher-\(h^{1,2}\) example: non-SUSY workflow at \(h^{1,2}=3\)

As in the SUSY workflow, the higher-dimensional non-SUSY check is useful for development but too expensive for a default tutorial pass. Set RUN_H12_3_EXAMPLES = True to run it locally.

if RUN_H12_3_EXAMPLES:
    # h12=3 example (model_ID=4 has hyperplanes + solutions near the sampling region).
    model_3_nonSUSY = jvc.FluxVacuaFinder(
        h12=3, model_ID=4, maximum_degree=2, map_to_fd=True)

    sampler_3_nonSUSY = jvc.data_sampler(
        model_3_nonSUSY,
        moduli_bounds=(2., 8.),
        dilaton_bounds=(np.sqrt(3) / 2, 10.),
        axion_bounds=(-0.5, 0.5),
        seed=42,
    )

    finder_3 = model_3_nonSUSY

    results_3 = finder_3.sample_critical_points(
        Q=50, noscale=True, moduli_max=500, sampler=sampler_3_nonSUSY, n_target=10,
        n_batch=100,
        max_batches=2,
        isd_mode="ISD-",
        solver="newton",
        step_size=0.5,
        newton_max_iters=300,
        verbose=True,
    )

    n3 = len(results_3)
    n3_ns = sum(1 for r in results_3 if not r['is_susy'])
    n3_min = sum(1 for r in results_3 if r['is_minimum'])
    print(f"\nh12=3 non-SUSY results: {n3} critical points "
          f"({n3_ns} non-SUSY, {n3_min} minima)")
else:
    results_3 = []
    print("RUN_H12_3_EXAMPLES = False — skipping the h12=3 non-SUSY example.")

Workflow D: ellipsoid-constrained tadpole-feasible seeds

Workflows A–C draw the input-half flux from a box and lean on the ISD solve and the tadpole constraint to discard infeasible configurations. At fixed \((z,\tau)\) the signed D3 tadpole \(f^{\mathsf T}\Sigma h = \tfrac12\, u^{\mathsf T} M(z,\tau)\, u\) is a positive-definite quadratic form in the input half \(u\), so the admissible set \(\{0 \le f^{\mathsf T}\Sigma h \le Q\}\) is an ellipsoid. Drawing \(u\) inside it — with tadpole_quadratic and sample_in_ellipsoid from jaxvacua.sampling (introduced in NB04) — yields tadpole-feasible seeds by construction. This is increasingly valuable at larger \(h^{1,2}\), where a box almost never lands in the admissible region. The feasible seeds then feed the same Newton refinement (find_solution_steps) used above.

# Workflow D: build tadpole-feasible flux seeds inside the ellipsoid, then refine.
from jaxvacua.sampling import tadpole_quadratic, sample_in_ellipsoid

Q, isd_mode, N = model.Q(), "F", 64

# Cone-interior, Kaehler-metric-positive moduli; one ellipsoid-feasible flux each.
zc, tauc, _ = sampler.initial_guesses(N, moduli_sampling_mode="cone",
                                      rns_key=jax.random.PRNGKey(0))
km = np.asarray(sampler.filter_by_km(zc, tauc))
zc, tauc = zc[km], tauc[km]
M = jax.vmap(lambda z, t: tadpole_quadratic(model, z, t, isd_mode))(zc, tauc)
keys = jax.random.split(jax.random.PRNGKey(1), zc.shape[0])
u = jax.vmap(lambda m, k: sample_in_ellipsoid(k, m, 0.7 * Q, 1))(M, keys)[:, 0, :]  # 0.7Q sub-bound
fluxes = sampler.ISD_sampling(zc, jnp.conj(zc), tauc, jnp.conj(tauc),
                              u.astype(jnp.complex128), mode=isd_mode,
                              output="full", return_integer_flux=True, vmap=True).real

# Keep tadpole-feasible, non-degenerate seeds (0 <= f^T Sigma h <= Q).
tad = np.asarray(model.tadpole(fluxes)).real
feas = (tad >= 0) & (tad <= Q) & np.any(np.asarray(fluxes) != 0, axis=1)
print(f"ellipsoid seeds: {100 * feas.mean():.0f}% tadpole-feasible ({int(feas.sum())}/{len(feas)})")

# Newton-refine the feasible seeds and count converged SUSY ISD vacua.
zc, tauc, fluxes = zc[feas], tauc[feas], fluxes[feas].astype(jnp.float64)
m1, t1, f1, _ = find_solution_steps(zc, tauc, fluxes)
for _ in range(12):
    m1, t1, f1, _ = find_solution_steps(m1, t1, f1)
res = np.sum(np.abs(np.asarray(
    model.DW(m1, jnp.conj(m1), t1, jnp.conj(t1), f1))), axis=1)
conv = np.isfinite(res) & (res < 1e-9)
print(f"converged ISD vacua (|DW| < 1e-9): {int(conv.sum())}/{len(res)}")

Take-aways

• NB07 is now the practical sampling notebook. It shows how to move from ISD-completed flux seeds to refined vacua without mixing in bounded enumeration or benchmark machinery.

• Workflow A exposes the moving parts. Use the manually wired, doubly-vmapped pipeline when you need to inspect intermediate ISD shifts or add custom logic between rounding and Newton refinement.

• Workflow B is the routine SUSY API. sample_SUSY_flux_vacua wraps the same pipeline behind a shorter call and should be the default for ordinary ensemble generation.

• Workflow C targets scalar-potential critical points. sample_critical_points solves \(\partial_I V = 0\) and returns classification metadata such as is_susy, is_minimum, V, |DW| and Nflux.

• Workflow D seeds the search inside the tadpole ellipsoid. tadpole_quadratic + sample_in_ellipsoid give input-half fluxes that satisfy \(0\le f^{\mathsf T}\Sigma h\le Q\) by construction, a scalable alternative to box sampling at larger \(h^{1,2}\).

• Diagnostics should appear immediately after sampling. Small moduli, \(W_0\), residual and tadpole plots catch bad regions or runaway-heavy runs before investing in larger scans.

Further reading

• NB04 — Sampling module — data_sampler infrastructure used by every workflow.

• NB05 — Finding flux vacua — Newton refinement and single-vacuum diagnostics.

• NB06 — ISD sampling principle — algorithmic intro to ISD completion.

• NB08 — Flux bounding — systematic bounded-flux enumeration and stochastic bounded sampling.

• NB16 — Sampling benchmarks — solver, sampler and tuning comparisons.

• NB17 — Flux-bounding case studies — literature-data validation against mirror-octic and Dataset B vacua.

• NB13 — Visualisation cookbook — publication-ready figures from sampling output.