diff --git a/.codecov.yml b/.codecov.yml
index 0e9aa6870b..17e0b61338 100644
--- a/.codecov.yml
+++ b/.codecov.yml
@@ -16,8 +16,11 @@ coverage:
     changes: false
     project:
       pymor:
+        # since we have no policy on this might as well turn this off entirely
+        target: 0%
         flags:
           - gitlab_ci
           - github_actions
     patch:
       pymor:
+        target: 0%
diff --git a/docs/source/substitutions.py b/docs/source/substitutions.py
index f7f9c2c436..c482fde83d 100644
--- a/docs/source/substitutions.py
+++ b/docs/source/substitutions.py
@@ -164,6 +164,8 @@
 
 .. |SymplecticBasis| replace:: :class:`~pymor.algorithms.symplectic.SymplecticBasis`
 .. |CanonicalSymplecticFormOperator| replace:: :class:`~pymor.operators.symplectic.CanonicalSymplecticFormOperator`
+
+.. |RNG| replace:: :class:`random number generator <pymor.tools.random.RNG>`
 '''
 
 substitutions = interfaces + common
diff --git a/src/pymor/algorithms/eigs.py b/src/pymor/algorithms/eigs.py
index 25609d88db..ae24ebde61 100644
--- a/src/pymor/algorithms/eigs.py
+++ b/src/pymor/algorithms/eigs.py
@@ -9,10 +9,11 @@
 from pymor.core.logger import getLogger
 from pymor.operators.constructions import IdentityOperator, InverseOperator
 from pymor.operators.interface import Operator
+from pymor.tools.random import new_rng
 
 
 def eigs(A, E=None, k=3, sigma=None, which='LM', b=None, l=None, maxiter=1000, tol=1e-13,
-         imag_tol=1e-12, complex_pair_tol=1e-12, complex_evp=False, left_evp=False, seed=0):
+         imag_tol=1e-12, complex_pair_tol=1e-12, complex_evp=False, left_evp=False):
     """Approximate a few eigenvalues of a linear |Operator|.
 
     Computes `k` eigenvalues `w` with corresponding eigenvectors `v` which solve
@@ -67,9 +68,6 @@ def eigs(A, E=None, k=3, sigma=None, which='LM', b=None, l=None, maxiter=1000, t
         are real setting this argument to `False` will increase stability and performance.
     left_evp
         If set to `True` compute left eigenvectors else compute right eigenvectors.
-    seed
-        Random seed which is used for computing the initial vector for the Arnoldi
-        iteration.
 
     Returns
     -------
@@ -93,7 +91,8 @@ def eigs(A, E=None, k=3, sigma=None, which='LM', b=None, l=None, maxiter=1000, t
         assert E.source == A.source
 
     if b is None:
-        b = A.source.random(seed=seed)
+        with new_rng(0):
+            b = A.source.random()
     else:
         assert b in A.source
 
diff --git a/src/pymor/algorithms/hapod.py b/src/pymor/algorithms/hapod.py
index a6500c610d..3eb7e0a752 100644
--- a/src/pymor/algorithms/hapod.py
+++ b/src/pymor/algorithms/hapod.py
@@ -10,6 +10,7 @@
 
 from pymor.algorithms.pod import pod
 from pymor.core.logger import getLogger
+from pymor.tools.random import spawn_rng
 
 
 class Node:
@@ -180,7 +181,7 @@ async def hapod_step(node):
 
         if node.children:
             modes, svals, snap_counts = zip(
-                *await asyncio.gather(*(hapod_step(c) for c in node.children))
+                *await asyncio.gather(*(spawn_rng(hapod_step(c)) for c in node.children))
             )
             for m, sv in zip(modes, svals):
                 m.scal(sv)
@@ -221,7 +222,7 @@ def main():
         nonlocal result
         result = asyncio.run(hapod_step(tree))
     result = None
-    hapod_thread = Thread(target=main)
+    hapod_thread = Thread(target=spawn_rng(main))
     hapod_thread.start()
     hapod_thread.join()
     return result
@@ -433,7 +434,7 @@ def __init__(self, executor, max_workers=None):
 
     def submit(self, f, *args):
         future = asyncio.get_event_loop().create_future()
-        self.queue.append((future, f, args))
+        self.queue.append((future, spawn_rng(f), args))
         asyncio.get_event_loop().create_task(self.run_task())
         return future
 
diff --git a/src/pymor/algorithms/lradi.py b/src/pymor/algorithms/lradi.py
index 878f55b0d3..1fecf52702 100644
--- a/src/pymor/algorithms/lradi.py
+++ b/src/pymor/algorithms/lradi.py
@@ -12,18 +12,17 @@
 from pymor.core.defaults import defaults
 from pymor.core.logger import getLogger
 from pymor.operators.constructions import IdentityOperator, InverseOperator
-from pymor.tools.random import get_random_state
+from pymor.tools.random import new_rng
 from pymor.vectorarrays.constructions import cat_arrays
 
 
 @defaults('lradi_tol', 'lradi_maxiter', 'lradi_shifts', 'projection_shifts_init_maxiter',
-          'projection_shifts_init_seed', 'projection_shifts_subspace_columns',
+          'projection_shifts_subspace_columns',
           'wachspress_large_ritz_num', 'wachspress_small_ritz_num', 'wachspress_tol')
 def lyap_lrcf_solver_options(lradi_tol=1e-10,
                              lradi_maxiter=500,
                              lradi_shifts='projection_shifts',
                              projection_shifts_init_maxiter=20,
-                             projection_shifts_init_seed=None,
                              projection_shifts_subspace_columns=6,
                              wachspress_large_ritz_num=50,
                              wachspress_small_ritz_num=25,
@@ -40,8 +39,6 @@ def lyap_lrcf_solver_options(lradi_tol=1e-10,
         See :func:`solve_lyap_lrcf`.
     projection_shifts_init_maxiter
         See :func:`projection_shifts_init`.
-    projection_shifts_init_seed
-        See :func:`projection_shifts_init`.
     projection_shifts_subspace_columns
         See :func:`projection_shifts`.
     wachspress_large_ritz_num
@@ -62,7 +59,6 @@ def lyap_lrcf_solver_options(lradi_tol=1e-10,
                       'shift_options':
                       {'projection_shifts': {'type': 'projection_shifts',
                                              'init_maxiter': projection_shifts_init_maxiter,
-                                             'init_seed': projection_shifts_init_seed,
                                              'subspace_columns': projection_shifts_subspace_columns},
                        'wachspress_shifts': {'type': 'wachspress_shifts',
                                              'large_ritz_num': wachspress_large_ritz_num,
@@ -194,14 +190,15 @@ def projection_shifts_init(A, E, B, shift_options):
     shifts
         A |NumPy array| containing a set of stable shift parameters.
     """
-    random_state = get_random_state(seed=shift_options['init_seed'])
+    rng = new_rng(0)
     for i in range(shift_options['init_maxiter']):
         Q = gram_schmidt(B, atol=0, rtol=0)
         shifts = spla.eigvals(A.apply2(Q, Q), E.apply2(Q, Q))
         shifts = shifts[shifts.real < 0]
         if shifts.size == 0:
             # use random subspace instead of span{B} (with same dimensions)
-            B = B.random(len(B), distribution='normal', random_state=random_state)
+            with rng:
+                B = B.random(len(B), distribution='normal')
         else:
             return shifts
     raise RuntimeError('Could not generate initial shifts for low-rank ADI iteration.')
diff --git a/src/pymor/algorithms/lrradi.py b/src/pymor/algorithms/lrradi.py
index 26d62c65d1..30917bce45 100644
--- a/src/pymor/algorithms/lrradi.py
+++ b/src/pymor/algorithms/lrradi.py
@@ -12,16 +12,15 @@
 from pymor.core.logger import getLogger
 from pymor.operators.constructions import IdentityOperator
 from pymor.algorithms.gram_schmidt import gram_schmidt
-from pymor.tools.random import get_random_state
+from pymor.tools.random import new_rng
 
 
 @defaults('lrradi_tol', 'lrradi_maxiter', 'lrradi_shifts', 'hamiltonian_shifts_init_maxiter',
-          'hamiltonian_shifts_init_seed', 'hamiltonian_shifts_subspace_columns')
+          'hamiltonian_shifts_subspace_columns')
 def ricc_lrcf_solver_options(lrradi_tol=1e-10,
                              lrradi_maxiter=500,
                              lrradi_shifts='hamiltonian_shifts',
                              hamiltonian_shifts_init_maxiter=20,
-                             hamiltonian_shifts_init_seed=None,
                              hamiltonian_shifts_subspace_columns=6):
     """Returns available Riccati equation solvers with default solver options.
 
@@ -35,8 +34,6 @@ def ricc_lrcf_solver_options(lrradi_tol=1e-10,
         See :func:`solve_ricc_lrcf`.
     hamiltonian_shifts_init_maxiter
         See :func:`hamiltonian_shifts_init`.
-    hamiltonian_shifts_init_seed
-        See :func:`hamiltonian_shifts_init`.
     hamiltonian_shifts_subspace_columns
         See :func:`hamiltonian_shifts`.
 
@@ -51,7 +48,6 @@ def ricc_lrcf_solver_options(lrradi_tol=1e-10,
                        'shift_options':
                        {'hamiltonian_shifts': {'type': 'hamiltonian_shifts',
                                                'init_maxiter': hamiltonian_shifts_init_maxiter,
-                                               'init_seed': hamiltonian_shifts_init_seed,
                                                'subspace_columns': hamiltonian_shifts_subspace_columns}}}}
 
 
@@ -229,7 +225,7 @@ def hamiltonian_shifts_init(A, E, B, C, shift_options):
     shifts
         A |NumPy array| containing a set of stable shift parameters.
     """
-    random_state = get_random_state(seed=shift_options['init_seed'])
+    rng = new_rng(0)
     for _ in range(shift_options['init_maxiter']):
         Q = gram_schmidt(C, atol=0, rtol=0)
         Ap = A.apply2(Q, Q)
@@ -249,7 +245,8 @@ def hamiltonian_shifts_init(A, E, B, C, shift_options):
         eigpairs = list(filter(lambda e: e[0].real < 0, eigpairs))
         if len(eigpairs) == 0:
             # use random subspace instead of span{C} (with same dimensions)
-            C = C.random(len(C), distribution='normal', random_state=random_state)
+            with rng:
+                C = C.random(len(C), distribution='normal')
             continue
         # find shift with most impact on convergence
         maxval = -1
diff --git a/src/pymor/algorithms/samdp.py b/src/pymor/algorithms/samdp.py
index d07d256bf2..7e3edff3f5 100644
--- a/src/pymor/algorithms/samdp.py
+++ b/src/pymor/algorithms/samdp.py
@@ -10,12 +10,13 @@
 from pymor.core.logger import getLogger
 from pymor.operators.constructions import IdentityOperator
 from pymor.operators.interface import Operator
+from pymor.tools.random import new_rng
 
 
 @defaults('which', 'tol', 'imagtol', 'conjtol', 'dorqitol', 'rqitol', 'maxrestart', 'krestart', 'init_shifts',
-          'rqi_maxiter', 'seed')
+          'rqi_maxiter')
 def samdp(A, E, B, C, nwanted, init_shifts=None, which='NR', tol=1e-10, imagtol=1e-6, conjtol=1e-8,
-          dorqitol=1e-4, rqitol=1e-10, maxrestart=100, krestart=20, rqi_maxiter=10, seed=0):
+          dorqitol=1e-4, rqitol=1e-10, maxrestart=100, krestart=20, rqi_maxiter=10):
     """Compute the dominant pole triplets and residues of the transfer function of an LTI system.
 
     This function uses the subspace accelerated dominant pole (SAMDP) algorithm as described in
@@ -70,8 +71,6 @@ def samdp(A, E, B, C, nwanted, init_shifts=None, which='NR', tol=1e-10, imagtol=
         Maximum dimension of search space before performing a restart.
     rqi_maxiter
         Maximum number of iterations for the two-sided Rayleigh quotient iteration.
-    seed
-        Random seed which is used for computing the initial shift and random restarts.
 
     Returns
     -------
@@ -106,7 +105,7 @@ def samdp(A, E, B, C, nwanted, init_shifts=None, which='NR', tol=1e-10, imagtol=
     k = 0
     nrestart = 0
     nr_converged = 0
-    np.random.seed(seed)
+    rng = new_rng(0)
 
     X = A.source.empty()
     Q = A.source.empty()
@@ -121,7 +120,7 @@ def samdp(A, E, B, C, nwanted, init_shifts=None, which='NR', tol=1e-10, imagtol=
     poles = np.empty(0)
 
     if init_shifts is None:
-        st = np.random.uniform() * 10.j
+        st = rng.uniform() * 10.j
         shift_nr = 0
         nr_shifts = 0
     else:
@@ -164,7 +163,7 @@ def samdp(A, E, B, C, nwanted, init_shifts=None, which='NR', tol=1e-10, imagtol=
         SH, UR, URt, res = _select_max_eig(H, G, X, V, B_defl, C_defl, which)
 
         if np.all(res < np.finfo(float).eps):
-            st = np.random.uniform() * 10.j
+            st = rng.uniform() * 10.j
             found = False
         else:
             found = True
@@ -287,7 +286,7 @@ def samdp(A, E, B, C, nwanted, init_shifts=None, which='NR', tol=1e-10, imagtol=
                     if found:
                         st = SH[0, 0]
                     else:
-                        st = np.random.uniform() * 10.j
+                        st = rng.uniform() * 10.j
 
                     if shift_nr < nr_shifts:
                         st = shifts[shift_nr]
diff --git a/src/pymor/basic.py b/src/pymor/basic.py
index ce05dd7625..f09f171db5 100644
--- a/src/pymor/basic.py
+++ b/src/pymor/basic.py
@@ -76,7 +76,7 @@
                                   SOBTReductor)
 from pymor.reductors.sor_irka import SORIRKAReductor
 
-from pymor.tools.random import default_random_state
+from pymor.tools.random import get_rng, new_rng
 
 from pymor.vectorarrays.constructions import cat_arrays
 from pymor.vectorarrays.list import ListVectorSpace
diff --git a/src/pymor/bindings/dunegdt.py b/src/pymor/bindings/dunegdt.py
index 2a53a5370f..e502dc5232 100644
--- a/src/pymor/bindings/dunegdt.py
+++ b/src/pymor/bindings/dunegdt.py
@@ -125,8 +125,8 @@ def real_zero_vector(self):
     def real_full_vector(self, value):
         return DuneXTVector(self.dune_vector_type(self.dim, value))
 
-    def real_random_vector(self, distribution, random_state, **kwargs):
-        values = _create_random_values(self.dim, distribution, random_state, **kwargs)
+    def real_random_vector(self, distribution, **kwargs):
+        values = _create_random_values(self.dim, distribution, **kwargs)
         return self.real_vector_from_numpy(values)
 
     def real_vector_from_numpy(self, data, ensure_copy=False):
diff --git a/src/pymor/bindings/fenics.py b/src/pymor/bindings/fenics.py
index 05e7f3b41f..7135a24502 100644
--- a/src/pymor/bindings/fenics.py
+++ b/src/pymor/bindings/fenics.py
@@ -157,9 +157,9 @@ def real_full_vector(self, value):
         impl += value
         return FenicsVector(impl)
 
-    def real_random_vector(self, distribution, random_state, **kwargs):
+    def real_random_vector(self, distribution, **kwargs):
         impl = df.Function(self.V).vector()
-        values = _create_random_values(impl.local_size(), distribution, random_state, **kwargs)
+        values = _create_random_values(impl.local_size(), distribution, **kwargs)
         impl[:] = np.ascontiguousarray(values)
         return FenicsVector(impl)
 
diff --git a/src/pymor/parallel/ipython.py b/src/pymor/parallel/ipython.py
index 765615f2a8..1564cf026c 100644
--- a/src/pymor/parallel/ipython.py
+++ b/src/pymor/parallel/ipython.py
@@ -147,7 +147,8 @@ def __init__(self, num_engines=None, **kwargs):
         else:
             self.view = self.client[:]
         self.logger.info(f'Connected to {len(self.view)} engines')
-        self.view.map_sync(_setup_worker, range(len(self.view)))
+        from pymor.tools.random import get_seed_seq
+        self.view.map_sync(_setup_worker, get_seed_seq().spawn(len(self.view)))
         self._remote_objects_created = Counter()
 
         if defaults.defaults_changes() > 0:
@@ -202,15 +203,13 @@ def _worker_call_function(function, loop, args, kwargs):
         return function(*args, **kwargs)
 
 
-def _setup_worker(worker_id):
+def _setup_worker(seed_seq):
     global _remote_objects
     _remote_objects = {}
-    # ensure that each worker starts with a different RandomState
-    from pymor.tools import random
-    import numpy as np
-    state = random.default_random_state()
-    new_state = np.random.RandomState(state.randint(0, 2**16) + worker_id)
-    random._default_random_state = new_state
+    # ensure that each worker starts with a different yet deterministically
+    # initialized rng
+    from pymor.tools.random import new_rng
+    new_rng(seed_seq).install()
 
 
 def _push_object(remote_id, obj):
diff --git a/src/pymor/parallel/mpi.py b/src/pymor/parallel/mpi.py
index fb6ff7cded..56150ba29f 100644
--- a/src/pymor/parallel/mpi.py
+++ b/src/pymor/parallel/mpi.py
@@ -8,6 +8,7 @@
 
 from pymor.parallel.basic import WorkerPoolBase
 from pymor.tools import mpi
+from pymor.tools.random import get_seed_seq
 
 
 class MPIPool(WorkerPoolBase):
@@ -18,6 +19,7 @@ def __init__(self):
         self.logger.info(f'Connected to {mpi.size} ranks')
         self._payload = mpi.call(mpi.function_call_manage, _setup_worker)
         self._apply(os.chdir, os.getcwd())
+        self._map(_setup_rng, [[[s] for s in get_seed_seq().spawn(mpi.size)]])
 
     def __del__(self):
         mpi.call(mpi.remove_object, self._payload)
@@ -90,15 +92,15 @@ def _worker_map_function(payload, function, **kwargs):
 
 
 def _setup_worker():
-    # ensure that each worker starts with a different RandomState
-    if not mpi.rank0:
-        from pymor.tools import random
-        import numpy as np
-        state = random.default_random_state()
-        new_state = np.random.RandomState(state.randint(0, 2**16) + mpi.rank)
-        random._default_random_state = new_state
     return [None]
 
 
+def _setup_rng(seed_seq):
+    # ensure that each worker starts with a different yet deterministically
+    # initialized rng
+    from pymor.tools.random import new_rng
+    new_rng(seed_seq).install()
+
+
 def _push_object(obj):
     return obj
diff --git a/src/pymor/parameters/base.py b/src/pymor/parameters/base.py
index 4dab1dd256..776629eb0a 100644
--- a/src/pymor/parameters/base.py
+++ b/src/pymor/parameters/base.py
@@ -21,7 +21,7 @@
 from pymor.tools.floatcmp import float_cmp_all
 from pymor.tools.frozendict import FrozenDict, SortedFrozenDict
 from pymor.tools.pprint import format_array
-from pymor.tools.random import get_random_state
+from pymor.tools.random import get_rng
 
 
 class Parameters(SortedFrozenDict):
@@ -554,7 +554,7 @@ def sample_uniformly(self, counts):
         return [Mu((k, np.array(v)) for k, v in zip(self.parameters, i))
                 for i in product(*iters)]
 
-    def sample_randomly(self, count=None, random_state=None, seed=None):
+    def sample_randomly(self, count=None):
         """Randomly sample |parameter values| from the space.
 
         Parameters
@@ -563,21 +563,12 @@ def sample_randomly(self, count=None, random_state=None, seed=None):
             If `None`, a single dict `mu` of |parameter values| is returned.
             Otherwise, the number of random samples to generate and return as
             a list of |parameter values| dicts.
-        random_state
-            :class:`~numpy.random.RandomState` to use for sampling.
-            If `None`, a new random state is generated using `seed`
-            as random seed, or the :func:`default <pymor.tools.random.default_random_state>`
-            random state is used.
-        seed
-            If not `None`, a new random state with this seed is used.
 
         Returns
         -------
         The sampled |parameter values|.
         """
-        assert not random_state or seed is None
-        random_state = get_random_state(random_state, seed)
-        get_param = lambda: Mu(((k, random_state.uniform(self.ranges[k][0], self.ranges[k][1], size))
+        get_param = lambda: Mu(((k, get_rng().uniform(self.ranges[k][0], self.ranges[k][1], size))
                                for k, size in self.parameters.items()))
         if count is None:
             return get_param()
diff --git a/src/pymor/reductors/h2.py b/src/pymor/reductors/h2.py
index 3f7acf5538..fbcc0636da 100644
--- a/src/pymor/reductors/h2.py
+++ b/src/pymor/reductors/h2.py
@@ -21,6 +21,7 @@
 from pymor.parameters.base import Mu
 from pymor.reductors.basic import LTIPGReductor
 from pymor.reductors.interpolation import LTIBHIReductor, TFBHIReductor
+from pymor.tools.random import new_rng
 
 
 class GenericIRKAReductor(BasicObject):
@@ -94,10 +95,10 @@ def _order_to_sigma_b_c(self, r):
         sigma = np.logspace(-1, 1, r)
         b = (np.ones((r, 1))
              if self.fom.dim_input == 1
-             else np.random.RandomState(0).normal(size=(r, self.fom.dim_input)))
+             else new_rng(0).normal(size=(r, self.fom.dim_input)))
         c = (np.ones((r, 1))
              if self.fom.dim_output == 1
-             else np.random.RandomState(0).normal(size=(r, self.fom.dim_output)))
+             else new_rng(0).normal(size=(r, self.fom.dim_output)))
         return sigma, b, c
 
     @staticmethod
diff --git a/src/pymor/reductors/neural_network.py b/src/pymor/reductors/neural_network.py
index 4fc3abb6a3..4af810b985 100644
--- a/src/pymor/reductors/neural_network.py
+++ b/src/pymor/reductors/neural_network.py
@@ -25,6 +25,7 @@
                                          NeuralNetworkStatefreeOutputModel,
                                          NeuralNetworkInstationaryModel,
                                          NeuralNetworkInstationaryStatefreeOutputModel)
+from pymor.tools.random import get_rng, get_seed_seq
 
 
 class NeuralNetworkReductor(BasicObject):
@@ -104,7 +105,7 @@ def reduce(self, hidden_layers='[(N+P)*3, (N+P)*3]', activation_function=torch.t
                loss_function=None, restarts=10, lr_scheduler=optim.lr_scheduler.StepLR,
                lr_scheduler_params={'step_size': 10, 'gamma': 0.7},
                es_scheduler_params={'patience': 10, 'delta': 0.}, weight_decay=0.,
-               log_loss_frequency=0, seed=0):
+               log_loss_frequency=0):
         """Reduce by training artificial neural networks.
 
         Parameters
@@ -154,9 +155,6 @@ def reduce(self, hidden_layers='[(N+P)*3, (N+P)*3]', activation_function=torch.t
             Frequency of epochs in which to log the current validation and
             training loss during training of the neural networks.
             If `0`, no intermediate logging of losses is done.
-        seed
-            Seed to use for various functions in PyTorch. Using a fixed seed,
-            it is possible to reproduce former results.
 
         Returns
         -------
@@ -169,9 +167,7 @@ def reduce(self, hidden_layers='[(N+P)*3, (N+P)*3]', activation_function=torch.t
         assert learning_rate > 0.
         assert weight_decay >= 0.
 
-        # set a seed for the PyTorch initialization of weights and biases
-        # and further PyTorch methods
-        torch.manual_seed(seed)
+        torch.manual_seed(get_seed_seq().spawn(1)[0].generate_state(1).item())
 
         # build a reduced basis using POD and compute training data
         if not hasattr(self, 'training_data'):
@@ -196,7 +192,7 @@ def reduce(self, hidden_layers='[(N+P)*3, (N+P)*3]', activation_function=torch.t
                 else:
                     number_validation_snapshots = int(len(self.training_data)*self.validation_ratio)
                     # randomly shuffle training data before splitting into two sets
-                    np.random.shuffle(self.training_data)
+                    get_rng().shuffle(self.training_data)
                     # split training data into validation and training set
                     self.validation_data = self.training_data[0:number_validation_snapshots]
                     self.training_data = self.training_data[number_validation_snapshots+1:]
@@ -231,7 +227,7 @@ def weighted_mse_loss_function(inputs, targets):
             self.neural_network, self.losses = multiple_restarts_training(self.training_data, self.validation_data,
                                                                           neural_network, target_loss, restarts,
                                                                           log_loss_frequency, training_parameters,
-                                                                          self.scaling_parameters, seed)
+                                                                          self.scaling_parameters)
 
         self._check_tolerances()
 
@@ -1082,7 +1078,7 @@ def closure():
 
 def multiple_restarts_training(training_data, validation_data, neural_network,
                                target_loss=None, max_restarts=10, log_loss_frequency=0,
-                               training_parameters={}, scaling_parameters={}, seed=None):
+                               training_parameters={}, scaling_parameters={}):
     """Algorithm that performs multiple restarts of neural network training.
 
     This method either performs a predefined number of restarts and returns
@@ -1114,9 +1110,6 @@ def multiple_restarts_training(training_data, validation_data, neural_network,
     scaling_parameters
         Additional parameters for scaling inputs respectively outputs,
         see :func:`train_neural_network` for more information.
-    seed
-        Seed to use for various functions in PyTorch. Using a fixed seed,
-        it is possible to reproduce former results.
 
     Returns
     -------
@@ -1136,10 +1129,7 @@ def multiple_restarts_training(training_data, validation_data, neural_network,
 
     logger = getLogger('pymor.algorithms.neural_network.multiple_restarts_training')
 
-    # if applicable, set a common seed for the PyTorch initialization
-    # of weights and biases and further PyTorch methods for all training runs
-    if seed:
-        torch.manual_seed(seed)
+    torch.manual_seed(get_seed_seq().spawn(1)[0].generate_state(1).item())
 
     # in case no training data is provided, return a neural network
     # that always returns zeros independent of the input
diff --git a/src/pymor/tools/random.py b/src/pymor/tools/random.py
index a07c02dc0e..88dd6f5921 100644
--- a/src/pymor/tools/random.py
+++ b/src/pymor/tools/random.py
@@ -2,54 +2,198 @@
 # Copyright pyMOR developers and contributors. All rights reserved.
 # License: BSD 2-Clause License (https://opensource.org/licenses/BSD-2-Clause)
 
+"""Methods for managing random state in pyMOR.
+
+Many algorithms potentially depend directly or indirectly on randomness.
+To ensure reproducible execution of pyMOR code without having to pass around
+a random number generator object everywhere, pyMOR manages a global |RNG|
+object. This object is initialized automatically from a configurable |default|
+random seed during startup and can be obtained by calling :func:`get_rng`.
+The returned object is a subclass of :class:`numpy.random.Generator` and
+inherits all its sampling methods.
+
+To locally reset the global |RNG| in order to deterministically sample random
+numbers independently of previously executed code, a new |RNG| can be created
+via :func:`new_rng` and installed by using it as a context manager. For
+instance, to sample a deterministic initialization vector for an iterative
+algorithm we can write:
+
+.. code:: python
+
+    with new_rng(12345):
+        U0 = some_operator.source.random()
+
+Using a single global random state can lead to either non-deterministic or
+correlated behavior in parallel or asynchronous code. :func:`get_rng` takes
+provisions to detect such situations and issue a warning. In such cases
+:func:`spawn_rng` needs to be called on the entry points of concurrent code
+paths to ensure the desired behavior. For an advanced example, see
+:mod:`pymor.algorithms.hapod`.
+"""
+
+from contextvars import ContextVar
+import inspect
+
 from pymor.core.defaults import defaults
 
 import numpy as np
 
 
-def get_random_state(random_state=None, seed=None):
-    """Returns a |NumPy| :class:`~numpy.random.RandomState`.
+def get_rng():
+    """Returns the current globally installed :class:`random number generator <RNG>`."""
+    rng_state = _get_rng_state()
+    rng_state[0] = True
+    _rng_state.set([False] + rng_state[1:])
+    return rng_state[1]
+
+
+@defaults('seed_seq')
+def new_rng(seed_seq=42):
+    """Creates a new |RNG| and returns it.
 
     Parameters
     ----------
-    random_state
-        If specified, this state is returned.
-    seed
-        If specified, the seed to initialize a new random state.
+    seed_seq
+        Entropy to seed the generator with. Either a :class:`~numpy.random.SeedSequence`
+        or an `int` or list of `ints` from which the :class:`~numpy.random.SeedSequence`
+        will be created. If `None`, entropy is sampled from the operating system.
 
     Returns
     -------
-    Either the provided, a newly created or the default `RandomState`
-    object.
+    The newly created |RNG|.
     """
-    assert random_state is None or seed is None
-    if random_state is not None:
-        return random_state
-    elif seed is not None:
-        return np.random.RandomState(seed)
-    else:
-        return default_random_state()
+    if not isinstance(seed_seq, np.random.SeedSequence):
+        seed_seq = np.random.SeedSequence(seed_seq)
+    return RNG(seed_seq)
+
+
+class RNG(np.random.Generator):
+    """Random number generator.
+
+    This class inherits from :class:`np.random.Generator` and inherits all its sampling
+    methods. Further, the class can be used as a context manager, which upon entry
+    installs the RNG as pyMOR's global RNG that is returned from :func:`get_rng`.
+    When the context is left, the previous global RNG is installed again.
+
+    When using a context manager is not feasible, i.e. in an interactive workflow, this
+    functionality can be accessed via the :meth:`~RNG.install` and :meth:`~RNG:uninstall`
+    methods.
+
+    A new instance of this class should be obtained using :func:`new_rng`.
+
+    Parameters
+    ----------
+    seed_seq
+        A :class:`~numpy.random.SeedSequence` to initialized the RNG with.
+    """
+
+    def __init__(self, seed_seq):
+        self.old_state = _rng_state.get(None)
+        super().__init__(np.random.default_rng(seed_seq).bit_generator)
+        self._seed_seq = seed_seq
 
+    def install(self):
+        """Installs the generator as pyMOR's global random generator."""
+        # Store a new rng together with seed_seq, as the latter cannot be recovered from the
+        # rng via a public interface (see https://github.com/numpy/numpy/issues/15322).
+        # The first field is a flag to indicate whether the current _rng_state has been consumed
+        # via get_rng. This is a safeguard to detect calls to get_rng in concurrent code
+        # paths.
+        _rng_state.set([False, self, self._seed_seq])
 
-@defaults('seed')
-def default_random_state(seed=42):
-    """Returns the default |NumPy| :class:`~numpy.random.RandomState`.
+    def uninstall(self):
+        """Restores the previously set global random generator."""
+        _rng_state.set(self.old_state)
+
+    def __enter__(self):
+        self.install()
+        return self
+
+    def __exit__(self, exc_type, exc_value, exc_tb):
+        self.uninstall()
+
+
+def get_seed_seq():
+    """Returns :class:`~np.random.SeedSequence` of the current global |RNG|.
+
+    This function returns the :class:`~np.random.SeedSequence` with which pyMOR's
+    currently installed global |RNG| has been initialized. The returned instance can
+    be used to deterministically create a new :class:`~np.random.SeedSequence` via
+    the :meth:`~np.random.SeedSequence.spawn` method, which then can be used to
+    initialize a new random generator in external library code or concurrent code
+    paths.
+    """
+    return _get_rng_state()[2]
+
+
+def spawn_rng(f):
+    """Wraps a function or coroutine to create a new |RNG| in concurrent code paths.
+
+    Calling this function on a function or coroutine object creates a wrapper which
+    will execute the wrapped function with a new globally installed |RNG|. This
+    ensures that random numbers in concurrent code paths (:mod:`threads <threading>`,
+    :mod:`multiprocessing`, :mod:`asyncio`) are deterministically generated yet
+    uncorrelated.
+
+    .. warning::
+        If the control flow within a single code path depends on communication events
+        with concurrent code, e.g., the order in which some parallel jobs finish,
+        deterministic behavior can no longer be guaranteed by just using :func:`spawn_rng`.
+        In such cases, the code additionally has to ensure that random numbers are sampled
+        independently of the communication order.
 
     Parameters
     ----------
-    seed
-        Seed to use for initializing the random state.
+    f
+        The function or coroutine to wrap.
 
     Returns
     -------
-    The default `RandomState` object.
+    The wrapped function or coroutine.
     """
-    global _default_random_state
+    seed_seq = get_seed_seq().spawn(1)[0]
+
+    if inspect.iscoroutine(f):
+
+        async def spawn_rng_wrapper():
+            with new_rng(seed_seq):
+                return await f
+
+        return spawn_rng_wrapper()
+
+    elif inspect.isfunction(f):
+
+        return _SpawnRngWrapper(f, seed_seq)  # use a class to obtain something picklable
+
+    else:
+        raise TypeError
+
+
+class _SpawnRngWrapper:
+    def __init__(self, f, seed_seq):
+        self.f, self.seed_seq = f, seed_seq
+
+    def __call__(self, *args, **kwargs):
+        with new_rng(self.seed_seq):
+            return self.f(*args, **kwargs)
 
-    if _default_random_state is None:
-        _default_random_state = np.random.RandomState(seed)
 
-    return _default_random_state
+def _get_rng_state():
+    try:
+        rng_state = _rng_state.get()
+    except LookupError:
+        import warnings
+        warnings.warn(
+            'get_rng called but _rng_state not initialized. (Call spawn_rng when creating new thread.) '
+            'Initializing a new RNG from the default seed. This may lead to correlated data.')
+        new_rng().install()
+        rng_state = _rng_state.get()
+    if rng_state[0]:
+        import warnings
+        warnings.warn('You are using the same RNG in concurrent code paths.\n'
+                      'This may lead to truly random, irreproducible behavior')
+    return rng_state
 
 
-_default_random_state = None
+_rng_state = ContextVar('_rng_state')
+new_rng().install()
diff --git a/src/pymor/tools/timing.py b/src/pymor/tools/timing.py
index 7a615a2e14..7592a39472 100644
--- a/src/pymor/tools/timing.py
+++ b/src/pymor/tools/timing.py
@@ -68,6 +68,7 @@ def new_func(*args, **kwargs):
 
 
 def busywait(amount):
+    from pymor.tools.random import get_rng
     arr = np.arange(1000)
     for _ in range(amount):
-        np.random.shuffle(arr)
+        get_rng().shuffle(arr)
diff --git a/src/pymor/vectorarrays/interface.py b/src/pymor/vectorarrays/interface.py
index 3bcaac83e3..10276c03ab 100644
--- a/src/pymor/vectorarrays/interface.py
+++ b/src/pymor/vectorarrays/interface.py
@@ -8,7 +8,7 @@
 
 from pymor.core.base import BasicObject, ImmutableObject, abstractmethod
 from pymor.core.defaults import defaults
-from pymor.tools.random import get_random_state
+from pymor.tools.random import get_rng
 
 
 class VectorArray(BasicObject):
@@ -146,11 +146,11 @@ def full(self, value, count=1, reserve=0):
         """
         return self.space.full(value, count, reserve=reserve)
 
-    def random(self, count=1, distribution='uniform', random_state=None, seed=None, reserve=0, **kwargs):
+    def random(self, count=1, distribution='uniform', reserve=0, **kwargs):
         """Create a |VectorArray| of vectors with random entries.
 
         This is a shorthand for
-        `self.space.random(count, distribution, radom_state, seed, **kwargs)`.
+        `self.space.random(count, distribution, **kwargs)`.
 
         Supported random distributions::
 
@@ -176,17 +176,10 @@ def random(self, count=1, distribution='uniform', random_state=None, seed=None,
             Mean for `'normal'` distribution (defaults to `0`).
         scale
             Standard deviation for `'normal'` distribution (defaults to `1`).
-        random_state
-            :class:`~numpy.random.RandomState` to use for sampling.
-            If `None`, a new random state is generated using `seed`
-            as random seed, or the :func:`default <pymor.tools.random.default_random_state>`
-            random state is used.
-        seed
-            If not `None`, a new radom state with this seed is used.
         reserve
             Hint for the backend to which length the array will grow.
         """
-        return self.space.random(count, distribution, random_state, seed, **kwargs)
+        return self.space.random(count, distribution, **kwargs)
 
     def empty(self, reserve=0):
         """Create an empty |VectorArray| of the same |VectorSpace|.
@@ -903,7 +896,7 @@ def full(self, value, count=1, reserve=0):
         """
         return self.from_numpy(np.full((count, self.dim), value))
 
-    def random(self, count=1, distribution='uniform', random_state=None, seed=None, reserve=0, **kwargs):
+    def random(self, count=1, distribution='uniform', reserve=0, **kwargs):
         """Create a |VectorArray| of vectors with random entries.
 
         Supported random distributions::
@@ -930,19 +923,10 @@ def random(self, count=1, distribution='uniform', random_state=None, seed=None,
             Mean for `'normal'` distribution (defaults to `0`).
         scale
             Standard deviation for `'normal'` distribution (defaults to `1`).
-        random_state
-            :class:`~numpy.random.RandomState` to use for sampling.
-            If `None`, a new random state is generated using `seed`
-            as random seed, or the :func:`default <pymor.tools.random.default_random_state>`
-            random state is used.
-        seed
-            If not `None`, a new random state with this seed is used.
         reserve
             Hint for the backend to which length the array will grow.
         """
-        assert random_state is None or seed is None
-        random_state = get_random_state(random_state, seed)
-        values = _create_random_values((count, self.dim), distribution, random_state, **kwargs)
+        values = _create_random_values((count, self.dim), distribution, **kwargs)
         return self.from_numpy(values)
 
     def empty(self, reserve=0):
@@ -995,10 +979,12 @@ def __hash__(self):
         return hash(self.id)
 
 
-def _create_random_values(shape, distribution, random_state, **kwargs):
+def _create_random_values(shape, distribution, **kwargs):
     if distribution not in ('uniform', 'normal'):
         raise NotImplementedError
 
+    rng = get_rng()
+
     if distribution == 'uniform':
         if not kwargs.keys() <= {'low', 'high'}:
             raise ValueError
@@ -1006,13 +992,13 @@ def _create_random_values(shape, distribution, random_state, **kwargs):
         high = kwargs.get('high', 1.)
         if high <= low:
             raise ValueError
-        return random_state.uniform(low, high, shape)
+        return rng.uniform(low, high, shape)
     elif distribution == 'normal':
         if not kwargs.keys() <= {'loc', 'scale'}:
             raise ValueError
         loc = kwargs.get('loc', 0.)
         scale = kwargs.get('scale', 1.)
-        return random_state.normal(loc, scale, shape)
+        return rng.normal(loc, scale, shape)
     else:
         assert False
 
diff --git a/src/pymor/vectorarrays/list.py b/src/pymor/vectorarrays/list.py
index a64d301aaa..6856fd1c69 100644
--- a/src/pymor/vectorarrays/list.py
+++ b/src/pymor/vectorarrays/list.py
@@ -5,7 +5,6 @@
 import numpy as np
 
 from pymor.core.base import BasicObject, abstractmethod, abstractclassmethod, classinstancemethod
-from pymor.tools.random import get_random_state
 from pymor.vectorarrays.interface import VectorArray, VectorArrayImpl, VectorSpace, _create_random_values
 
 
@@ -536,8 +535,8 @@ def ones_vector(self):
     def full_vector(self, value):
         return self.vector_from_numpy(np.full(self.dim, value))
 
-    def random_vector(self, distribution, random_state, **kwargs):
-        values = _create_random_values(self.dim, distribution, random_state, **kwargs)
+    def random_vector(self, distribution, **kwargs):
+        values = _create_random_values(self.dim, distribution, **kwargs)
         return self.vector_from_numpy(values)
 
     @abstractmethod
@@ -567,13 +566,11 @@ def full(self, value, count=1, reserve=0):
         assert count >= 0 and reserve >= 0
         return ListVectorArray(self, ListVectorArrayImpl([self.full_vector(value) for _ in range(count)], self))
 
-    def random(self, count=1, distribution='uniform', random_state=None, seed=None, reserve=0, **kwargs):
+    def random(self, count=1, distribution='uniform', reserve=0, **kwargs):
         assert count >= 0 and reserve >= 0
-        assert random_state is None or seed is None
-        random_state = get_random_state(random_state, seed)
         return ListVectorArray(
             self,
-            ListVectorArrayImpl([self.random_vector(distribution=distribution, random_state=random_state, **kwargs)
+            ListVectorArrayImpl([self.random_vector(distribution=distribution, **kwargs)
                                  for _ in range(count)], self)
         )
 
@@ -621,12 +618,12 @@ def real_full_vector(self, value):
     def full_vector(self, value):
         return self.vector_type(self.real_full_vector(value), None)
 
-    def real_random_vector(self, distribution, random_state, **kwargs):
-        values = _create_random_values(self.dim, distribution, random_state, **kwargs)
+    def real_random_vector(self, distribution, **kwargs):
+        values = _create_random_values(self.dim, distribution, **kwargs)
         return self.real_vector_from_numpy(values)
 
-    def random_vector(self, distribution, random_state, **kwargs):
-        return self.vector_type(self.real_random_vector(distribution, random_state, **kwargs), None)
+    def random_vector(self, distribution, **kwargs):
+        return self.vector_type(self.real_random_vector(distribution, **kwargs), None)
 
     @abstractmethod
     def real_make_vector(self, obj):
diff --git a/src/pymor/vectorarrays/numpy.py b/src/pymor/vectorarrays/numpy.py
index 322181156f..a98c0fa351 100644
--- a/src/pymor/vectorarrays/numpy.py
+++ b/src/pymor/vectorarrays/numpy.py
@@ -8,7 +8,6 @@
 from scipy.sparse import issparse
 
 from pymor.core.base import classinstancemethod
-from pymor.tools.random import get_random_state
 from pymor.vectorarrays.interface import VectorArray, VectorArrayImpl, VectorSpace, _create_random_values
 
 
@@ -237,13 +236,11 @@ def full(self, value, count=1, reserve=0):
         assert reserve >= 0
         return NumpyVectorArray(self, NumpyVectorArrayImpl(np.full((max(count, reserve), self.dim), value), count))
 
-    def random(self, count=1, distribution='uniform', random_state=None, seed=None, reserve=0, **kwargs):
+    def random(self, count=1, distribution='uniform', reserve=0, **kwargs):
         assert count >= 0
         assert reserve >= 0
-        assert random_state is None or seed is None
-        random_state = get_random_state(random_state, seed)
         va = self.zeros(count, reserve)
-        va.impl._array[:count] = _create_random_values((count, self.dim), distribution, random_state, **kwargs)
+        va.impl._array[:count] = _create_random_values((count, self.dim), distribution, **kwargs)
         return va
 
     @classinstancemethod
diff --git a/src/pymordemos/dmd_identification.py b/src/pymordemos/dmd_identification.py
index 88ca3177ce..2b5ff1a24a 100755
--- a/src/pymordemos/dmd_identification.py
+++ b/src/pymordemos/dmd_identification.py
@@ -13,12 +13,9 @@
 
 def main(
         n: int = Option(4, help='Dimension of the state.'),
-        M: int = Option(10, help='Number of data pairs.'),
-        seed: int = Option(42, help='Random seed.')
+        M: int = Option(10, help='Number of data pairs.')
 ):
-    np.random.seed(seed)
-
-    A = np.random.rand(n, n)
+    A = get_rng().random((n, n))
     A = A / np.linalg.norm(A)
     print(f'A: {A}')
     X = np.zeros((M + 1, n))
diff --git a/src/pymordemos/output_error_estimation.py b/src/pymordemos/output_error_estimation.py
index fb5bbf652c..9b53cb6625 100755
--- a/src/pymordemos/output_error_estimation.py
+++ b/src/pymordemos/output_error_estimation.py
@@ -37,9 +37,9 @@ def main(
         # an output which is actually a lincomb operator
         fom = create_fom(grid_intervals, vector_valued_output=True)
         dim_source = fom.output_functional.source.dim
-        np.random.seed(1)
-        random_matrix_1 = np.random.rand(2, dim_source)
-        random_matrix_2 = np.random.rand(2, dim_source)
+        rng = get_rng()
+        random_matrix_1 = rng.random((2, dim_source))
+        random_matrix_2 = rng.random((2, dim_source))
         op1 = NumpyMatrixOperator(random_matrix_1, source_id='STATE')
         op2 = NumpyMatrixOperator(random_matrix_2, source_id='STATE')
         ops = [op1, op2]
@@ -52,7 +52,7 @@ def main(
         if fom_number == 3:
             fom = fom.with_(output_functional=fom.rhs.operators[0].H)
         else:
-            random_matrix_1 = np.random.rand(2, fom.solution_space.dim)
+            random_matrix_1 = get_rng().random((2, fom.solution_space.dim))
             op = NumpyMatrixOperator(random_matrix_1, source_id='STATE')
             fom = fom.with_(output_functional=op)
 
diff --git a/src/pymordemos/parabolic_mor.py b/src/pymordemos/parabolic_mor.py
index e179d5ba0f..997318f44c 100755
--- a/src/pymordemos/parabolic_mor.py
+++ b/src/pymordemos/parabolic_mor.py
@@ -66,7 +66,7 @@ def main(
         rom, fom=fom, reductor=reductor, error_estimator=True,
         error_norms=[lambda U: DT * np.sqrt(np.sum(fom.h1_0_semi_norm(U)[1:]**2))],
         error_norm_names=['l^2-h^1'],
-        condition=False, test_mus=parameter_space.sample_randomly(test, seed=999)
+        condition=False, test_mus=parameter_space.sample_randomly(test)
     )
 
     # show results
diff --git a/src/pymordemos/thermalblock.py b/src/pymordemos/thermalblock.py
index a09c17f455..721d703a0b 100755
--- a/src/pymordemos/thermalblock.py
+++ b/src/pymordemos/thermalblock.py
@@ -169,7 +169,7 @@ def main(
                                        error_estimator=True,
                                        error_norms=(fom.h1_0_semi_norm, fom.l2_norm),
                                        condition=True,
-                                       test_mus=parameter_space.sample_randomly(test, seed=999),
+                                       test_mus=parameter_space.sample_randomly(test),
                                        basis_sizes=0 if plot_error_sequence else 1,
                                        pool=None if fenics else pool  # cannot pickle FEniCS model
                                        )
diff --git a/src/pymortests/algorithms/ei.py b/src/pymortests/algorithms/ei.py
index 47507948dc..7f34aa61cd 100644
--- a/src/pymortests/algorithms/ei.py
+++ b/src/pymortests/algorithms/ei.py
@@ -7,6 +7,7 @@
 from pymor.algorithms.pod import pod
 from pymor.operators.ei import EmpiricalInterpolatedOperator
 from pymor.reductors.basic import StationaryRBReductor
+from pymor.tools.random import new_rng
 from pymortests.base import runmodule, assert_all_almost_equal
 
 
@@ -18,7 +19,9 @@ def test_ei_restricted_to_full(stationary_models):
     ei_op = EmpiricalInterpolatedOperator(op, collateral_basis=cb, interpolation_dofs=dofs, triangular=True)
     ei_model = model.with_(operator=ei_op)
 
-    for mu in model.parameters.space(1, 2).sample_randomly(3, seed=234):
+    with new_rng(234):
+        mus = model.parameters.space(1, 2).sample_randomly(3)
+    for mu in mus:
         a = model.solve(mu)
         b = ei_model.solve(mu)
         assert_all_almost_equal(a, b, rtol=1e-4)
@@ -41,7 +44,9 @@ def test_ei_rom(stationary_models):
 
     U = fom.solution_space.empty()
     base_mus = []
-    for mu in fom.parameters.space(1, 2).sample_randomly(3, seed=234):
+    with new_rng(234):
+        mus = fom.parameters.space(1, 2).sample_randomly(3)
+    for mu in mus:
         a = fom.solve(mu)
         U.append(a)
         base_mus.append(mu)
diff --git a/src/pymortests/algorithms/symplectic.py b/src/pymortests/algorithms/symplectic.py
index ba0c804370..674a7e4288 100644
--- a/src/pymortests/algorithms/symplectic.py
+++ b/src/pymortests/algorithms/symplectic.py
@@ -5,6 +5,7 @@
                                          psd_svd_like_decomp,
                                          symplectic_gram_schmidt)
 from pymor.operators.symplectic import CanonicalSymplecticFormOperator
+from pymor.tools.random import new_rng
 from pymor.vectorarrays.block import BlockVectorSpace
 from pymor.vectorarrays.numpy import NumpyVectorSpace
 
@@ -23,7 +24,8 @@ def test_symplecticity(key_method):
     half_space = NumpyVectorSpace(half_dim)
     phase_space = BlockVectorSpace([half_space] * 2)
     n_data = 100
-    U = phase_space.random(n_data, seed=42)
+    with new_rng(42):
+        U = phase_space.random(n_data)
     modes = 10
     basis = METHODS_DICT[key_method](U, modes)
     tis_basis = basis.transposed_symplectic_inverse()
@@ -37,7 +39,8 @@ def test_orthonormality(key_orthosympl_method):
     half_space = NumpyVectorSpace(half_dim)
     phase_space = BlockVectorSpace([half_space] * 2)
     n_data = 100
-    U = phase_space.random(n_data, seed=42)
+    with new_rng(42):
+        U = phase_space.random(n_data)
     modes = 10
     basis = METHODS_DICT[key_orthosympl_method](U, modes).to_array()
     assert np.allclose(basis.inner(basis), np.eye(modes))
@@ -53,13 +56,15 @@ def test_symplectic_gram_schmidt(test_orthonormality, reiterate):
     J = CanonicalSymplecticFormOperator(phase_space)
     half_red_dim = 10
 
-    E = phase_space.random(half_red_dim, seed=42)
+    with new_rng(42):
+        E = phase_space.random(half_red_dim)
     if test_orthonormality:
         # special choice, such that result is orthosymplectic
         F = J.apply(E)
     else:
         # less structure in snapshots, no orthogonality
-        F = phase_space.random(half_red_dim, seed=43)
+        with new_rng(43):
+            F = phase_space.random(half_red_dim)
 
     S, Lambda = symplectic_gram_schmidt(E, F, return_Lambda=True, reiterate=reiterate)
     tsi_S = S.transposed_symplectic_inverse()
diff --git a/src/pymortests/demos.py b/src/pymortests/demos.py
index 342d293718..01a3647e22 100644
--- a/src/pymortests/demos.py
+++ b/src/pymortests/demos.py
@@ -230,7 +230,7 @@ def nop(*args, **kwargs):
 
     # reset default RandomState
     import pymor.tools.random
-    pymor.tools.random._default_random_state = None
+    pymor.tools.random.new_rng().install()
 
     result = None
     try:
diff --git a/src/pymortests/fixtures/operator.py b/src/pymortests/fixtures/operator.py
index e26f43e8a2..c8da149088 100644
--- a/src/pymortests/fixtures/operator.py
+++ b/src/pymortests/fixtures/operator.py
@@ -11,6 +11,7 @@
 from pymor.operators.interface import Operator
 from pymor.operators.list import NumpyListVectorArrayMatrixOperator
 from pymor.operators.numpy import NumpyMatrixOperator
+from pymor.tools.random import new_rng
 from pymor.vectorarrays.numpy import NumpyVectorSpace
 
 
@@ -137,7 +138,9 @@ def thermalblock_factory(xblocks, yblocks, diameter, seed):
         U.append(iop.as_vector(f.parameters.parse(exp)))
     for exp in np.random.random(6):
         V.append(iop.as_vector(f.parameters.parse(exp)))
-    return m.operator, p.parameter_space.sample_randomly(1, seed=seed)[0], U, V, m.h1_product, m.l2_product
+    with new_rng(seed):
+        mu = p.parameter_space.sample_randomly(1)[0]
+    return m.operator, mu, U, V, m.h1_product, m.l2_product
 
 
 def thermalblock_assemble_factory(xblocks, yblocks, diameter, seed):
diff --git a/src/pymortests/low_rank_op.py b/src/pymortests/low_rank_op.py
index b53c050479..d61d798335 100644
--- a/src/pymortests/low_rank_op.py
+++ b/src/pymortests/low_rank_op.py
@@ -8,19 +8,20 @@
 from pymor.algorithms.lincomb import assemble_lincomb
 from pymor.operators.constructions import LowRankOperator, LowRankUpdatedOperator
 from pymor.operators.numpy import NumpyMatrixOperator
+from pymor.tools.random import new_rng
 from pymor.vectorarrays.numpy import NumpyVectorSpace
 
 
 def construct_operators_and_vectorarrays(m, n, r, k, seed=0):
     space_m = NumpyVectorSpace(m)
     space_n = NumpyVectorSpace(n)
-    rng = np.random.RandomState(seed)
-    A = NumpyMatrixOperator(rng.randn(m, n))
-    L = space_m.random(r, distribution='normal', random_state=rng)
-    C = rng.randn(r, r)
-    R = space_n.random(r, distribution='normal', random_state=rng)
-    U = space_n.random(k, distribution='normal', random_state=rng)
-    V = space_m.random(k, distribution='normal', random_state=rng)
+    with new_rng(seed) as rng:
+        A = NumpyMatrixOperator(rng.normal(size=(m, n)))
+        L = space_m.random(r, distribution='normal')
+        C = rng.normal(size=(r, r))
+        R = space_n.random(r, distribution='normal')
+        U = space_n.random(k, distribution='normal')
+        V = space_m.random(k, distribution='normal')
     return A, L, C, R, U, V
 
 
diff --git a/src/pymortests/model.py b/src/pymortests/model.py
index 3243d9cba8..1d88729594 100644
--- a/src/pymortests/model.py
+++ b/src/pymortests/model.py
@@ -14,6 +14,7 @@
 from pymor.models.symplectic import QuadraticHamiltonianModel
 from pymor.operators.block import BlockDiagonalOperator
 from pymor.operators.constructions import IdentityOperator
+from pymor.tools.random import new_rng
 from pymor.vectorarrays.numpy import NumpyVectorSpace
 from pymortests.base import runmodule
 from pymortests.pickling import assert_picklable, assert_picklable_without_dumps_function
@@ -32,7 +33,9 @@ def test_pickle_by_solving(model):
     m2 = loads(dumps(m))
     m.disable_caching()
     m2.disable_caching()
-    for mu in m.parameters.space(1, 2).sample_randomly(3, seed=234):
+    with new_rng(234):
+        mus = m.parameters.space(1, 2).sample_randomly(3)
+    for mu in mus:
         assert np.all(almost_equal(m.solve(mu), m2.solve(mu)))
 
 
diff --git a/src/pymortests/rand_la.py b/src/pymortests/rand_la.py
index 881362a62b..283c6bd72a 100644
--- a/src/pymortests/rand_la.py
+++ b/src/pymortests/rand_la.py
@@ -9,6 +9,7 @@
 from pymor.algorithms.rand_la import rrf, adaptive_rrf, random_ghep, random_generalized_svd
 from pymor.operators.numpy import NumpyMatrixOperator
 from pymor.operators.constructions import VectorArrayOperator
+from pymor.tools.random import new_rng
 
 
 def test_adaptive_rrf():
@@ -22,10 +23,14 @@ def test_adaptive_rrf():
     B = B.dot(B.T)
     source_product = NumpyMatrixOperator(B)
 
-    C = range_product.range.random(10, seed=10)
+    with new_rng(10):
+        C = range_product.range.random(10)
     op = VectorArrayOperator(C)
 
-    D = range_product.range.random(10, seed=11)+1j*range_product.range.random(10, seed=12)
+    with new_rng(11):
+        D = range_product.range.random(10)
+    with new_rng(12):
+        D += 1j*range_product.range.random(10)
     op_complex = VectorArrayOperator(D)
 
     Q1 = adaptive_rrf(op, source_product, range_product)
@@ -47,10 +52,14 @@ def test_rrf():
     B = B @ B.T
     source_product = NumpyMatrixOperator(B)
 
-    C = range_product.range.random(10, seed=10)
+    with new_rng(10):
+        C = range_product.range.random(10)
     op = VectorArrayOperator(C)
 
-    D = range_product.range.random(10, seed=11)+1j*range_product.range.random(10, seed=12)
+    with new_rng(11):
+        D = range_product.range.random(10)
+    with new_rng(12):
+        D += 1j*range_product.range.random(10)
     op_complex = VectorArrayOperator(D)
 
     Q1 = rrf(op, source_product, range_product)
diff --git a/src/pymortests/testdata/check_results/test_parabolic_mor_results/ca3b1109e5bd10d21b0b9bc432d6ac101cf1fa4d b/src/pymortests/testdata/check_results/test_parabolic_mor_results/ca3b1109e5bd10d21b0b9bc432d6ac101cf1fa4d
index 4583a8dbd5..01a375f087 100644
Binary files a/src/pymortests/testdata/check_results/test_parabolic_mor_results/ca3b1109e5bd10d21b0b9bc432d6ac101cf1fa4d and b/src/pymortests/testdata/check_results/test_parabolic_mor_results/ca3b1109e5bd10d21b0b9bc432d6ac101cf1fa4d differ
diff --git a/src/pymortests/testdata/check_results/test_thermalblock_results/105686ef68a9b4b3655c790d395526c623efd162 b/src/pymortests/testdata/check_results/test_thermalblock_results/105686ef68a9b4b3655c790d395526c623efd162
index 2e2c6a2583..02813ca032 100644
Binary files a/src/pymortests/testdata/check_results/test_thermalblock_results/105686ef68a9b4b3655c790d395526c623efd162 and b/src/pymortests/testdata/check_results/test_thermalblock_results/105686ef68a9b4b3655c790d395526c623efd162 differ
diff --git a/src/pymortests/testdata/check_results/test_thermalblock_results/29e5f263a02a275f4d29b7e5bdf32ad20c1e6346 b/src/pymortests/testdata/check_results/test_thermalblock_results/29e5f263a02a275f4d29b7e5bdf32ad20c1e6346
index 4c45447997..1a8cdb3a54 100644
Binary files a/src/pymortests/testdata/check_results/test_thermalblock_results/29e5f263a02a275f4d29b7e5bdf32ad20c1e6346 and b/src/pymortests/testdata/check_results/test_thermalblock_results/29e5f263a02a275f4d29b7e5bdf32ad20c1e6346 differ
diff --git a/src/pymortests/testdata/check_results/test_thermalblock_results/38a8bfe939551df7383a6058a2862fb9d2eea5a4 b/src/pymortests/testdata/check_results/test_thermalblock_results/38a8bfe939551df7383a6058a2862fb9d2eea5a4
index ee85571d15..652e715aff 100644
Binary files a/src/pymortests/testdata/check_results/test_thermalblock_results/38a8bfe939551df7383a6058a2862fb9d2eea5a4 and b/src/pymortests/testdata/check_results/test_thermalblock_results/38a8bfe939551df7383a6058a2862fb9d2eea5a4 differ
diff --git a/src/pymortests/testdata/check_results/test_thermalblock_results/a92ec4631a6f0edc2eae1cc3b0fa17af6c467cef b/src/pymortests/testdata/check_results/test_thermalblock_results/a92ec4631a6f0edc2eae1cc3b0fa17af6c467cef
index 65d0da2f97..4f9176cdc2 100644
Binary files a/src/pymortests/testdata/check_results/test_thermalblock_results/a92ec4631a6f0edc2eae1cc3b0fa17af6c467cef and b/src/pymortests/testdata/check_results/test_thermalblock_results/a92ec4631a6f0edc2eae1cc3b0fa17af6c467cef differ
diff --git a/src/pymortests/testdata/check_results/test_thermalblock_results/c8a460e412d9170f12508ab3608fd7e1e1750603 b/src/pymortests/testdata/check_results/test_thermalblock_results/c8a460e412d9170f12508ab3608fd7e1e1750603
index a54e1b621e..cd210e98d4 100644
Binary files a/src/pymortests/testdata/check_results/test_thermalblock_results/c8a460e412d9170f12508ab3608fd7e1e1750603 and b/src/pymortests/testdata/check_results/test_thermalblock_results/c8a460e412d9170f12508ab3608fd7e1e1750603 differ
diff --git a/src/pymortests/testdata/check_results/test_thermalblock_results/ce0770de067386276c9e89afa2d63daeb512a9f5 b/src/pymortests/testdata/check_results/test_thermalblock_results/ce0770de067386276c9e89afa2d63daeb512a9f5
index b4616c0f11..c6ac51d892 100644
Binary files a/src/pymortests/testdata/check_results/test_thermalblock_results/ce0770de067386276c9e89afa2d63daeb512a9f5 and b/src/pymortests/testdata/check_results/test_thermalblock_results/ce0770de067386276c9e89afa2d63daeb512a9f5 differ
diff --git a/src/pymortests/testdata/check_results/test_thermalblock_results/ea8f66e263715fd0185d0dd23ab6d57feb6ee898 b/src/pymortests/testdata/check_results/test_thermalblock_results/ea8f66e263715fd0185d0dd23ab6d57feb6ee898
index 91a76bd975..f4725c3c4a 100644
Binary files a/src/pymortests/testdata/check_results/test_thermalblock_results/ea8f66e263715fd0185d0dd23ab6d57feb6ee898 and b/src/pymortests/testdata/check_results/test_thermalblock_results/ea8f66e263715fd0185d0dd23ab6d57feb6ee898 differ
diff --git a/src/pymortests/vectorarray.py b/src/pymortests/vectorarray.py
index fcfb7119ee..4fbd467516 100644
--- a/src/pymortests/vectorarray.py
+++ b/src/pymortests/vectorarray.py
@@ -14,6 +14,7 @@
 from pymor.vectorarrays.interface import VectorSpace
 from pymor.vectorarrays.numpy import NumpyVectorSpace
 from pymor.tools.floatcmp import float_cmp, bounded
+from pymor.tools.random import new_rng
 from pymortests.base import might_exceed_deadline
 from pymortests.pickling import assert_picklable_without_dumps_function
 import pymortests.strategies as pyst
@@ -166,7 +167,8 @@ def _test_random_uniform(vector_array, realizations, low, high):
         return
     seed = 123
     try:
-        v = vector_array.random(c, low=low, high=high, seed=seed)
+        with new_rng(seed):
+            v = vector_array.random(c, low=low, high=high)
     except ValueError as e:
         if high <= low:
             return
@@ -182,7 +184,8 @@ def _test_random_uniform(vector_array, realizations, low, high):
         assert np.all(x >= low)
     except NotImplementedError:
         pass
-    vv = vector_array.random(c, distribution='uniform', low=low, high=high, seed=seed)
+    with new_rng(seed):
+        vv = vector_array.random(c, distribution='uniform', low=low, high=high)
     assert np.allclose((v - vv).sup_norm(), 0.)
 
 
@@ -199,7 +202,8 @@ def test_random_normal(vector_array, realizations, loc, scale):
         return
     seed = 123
     try:
-        v = vector_array.random(c, 'normal', loc=loc, scale=scale, seed=seed)
+        with new_rng(seed):
+            v = vector_array.random(c, 'normal', loc=loc, scale=scale)
     except ValueError as e:
         if scale <= 0:
             return
@@ -222,7 +226,8 @@ def test_random_normal(vector_array, realizations, loc, scale):
         bounded(lower, upper, loc)
     except NotImplementedError:
         pass
-    vv = vector_array.random(c, 'normal', loc=loc, scale=scale, seed=seed)
+    with new_rng(seed):
+        vv = vector_array.random(c, 'normal', loc=loc, scale=scale)
     data = vv.to_numpy()
     # due to scaling data might actually now include nan or inf
     assume(not np.isnan(data).any())