diff --git a/tests/test_active_plasma_lens.py b/tests/test_active_plasma_lens.py
index d86f794..be80e6a 100644
--- a/tests/test_active_plasma_lens.py
+++ b/tests/test_active_plasma_lens.py
@@ -29,7 +29,7 @@ def test_active_plasma_lens():
     apl.track(bunch)
     bunch_params = analyze_bunch(bunch)
     gamma_x = bunch_params['gamma_x']
-    assert approx(gamma_x, rel=1e-10) == 92.38646379897074
+    assert approx(gamma_x, rel=1e-10) == 92.38407675999406
 
 
 def test_active_plasma_lens_with_wakefields():
@@ -64,7 +64,7 @@ def test_active_plasma_lens_with_wakefields():
     # Analyze and check results.
     bunch_params = analyze_bunch(bunch)
     gamma_x = bunch_params['gamma_x']
-    assert approx(gamma_x, rel=1e-10) == 77.31995824746237
+    assert approx(gamma_x, rel=1e-10) == 77.32021188373825
 
 
 if __name__ == '__main__':
diff --git a/tests/test_beamline.py b/tests/test_beamline.py
index 62ca47a..44d0605 100644
--- a/tests/test_beamline.py
+++ b/tests/test_beamline.py
@@ -30,12 +30,12 @@ def test_single_element():
     bl = Beamline([deepcopy(plasma)])
 
     # Track plasma.
-    bunch_1 = deepcopy(bunch)
+    bunch_1 = bunch.copy()
     out_dir_1 = os.path.join(output_folder, 'plasma_diags')
     plasma.track(bunch_1, opmd_diag=True, diag_dir=out_dir_1)
 
     # Track beamline.
-    bunch_2 = deepcopy(bunch)
+    bunch_2 = bunch.copy()
     out_dir_2 = os.path.join(output_folder, 'bl_diags')
     bl.track(bunch_2, opmd_diag=True, diag_dir=out_dir_2)
 
@@ -81,7 +81,7 @@ def test_multiple_element():
     bl = Beamline([deepcopy(d1), deepcopy(plasma), deepcopy(d2)])
 
     # Track elements individually.
-    bunch_1 = deepcopy(bunch)
+    bunch_1 = bunch.copy()
     out_dir_d1 = os.path.join(output_folder, 'd1_diags')
     out_dir_plasma = os.path.join(output_folder, 'plasma_diags')
     out_dir_d2 = os.path.join(output_folder, 'd2_diags')
@@ -90,7 +90,7 @@ def test_multiple_element():
     d2.track(bunch_1, opmd_diag=True, diag_dir=out_dir_d2)
 
     # Track beamline.
-    bunch_2 = deepcopy(bunch)
+    bunch_2 = bunch.copy()
     out_dir_bl = os.path.join(output_folder, 'bl_diags')
     bl.track(bunch_2, opmd_diag=True, diag_dir=out_dir_bl)
 
diff --git a/tests/test_custom_blowout.py b/tests/test_custom_blowout.py
index d3d9792..6752986 100644
--- a/tests/test_custom_blowout.py
+++ b/tests/test_custom_blowout.py
@@ -46,7 +46,7 @@ def test_custom_blowout_wakefield(make_plots=False):
 
     bunch_params = analyze_bunch(bunch)
     rel_ene_sp = bunch_params['rel_ene_spread']
-    assert approx(rel_ene_sp, rel=1e-10) == 0.21192488237458038
+    assert approx(rel_ene_sp, rel=1e-10) == 0.21192494153185745
 
     if make_plots:
         # Analyze bunch evolution.
diff --git a/tests/test_field_element.py b/tests/test_field_element.py
index 1cd65f4..b1841dc 100644
--- a/tests/test_field_element.py
+++ b/tests/test_field_element.py
@@ -102,7 +102,7 @@ def test_field_element_error():
     with raises(ValueError) as e_info:
         element.track(bunch, opmd_diag=False)
     # This one should instead work.
-    element.track([bunch, copy.deepcopy(bunch)], opmd_diag=False)
+    element.track([bunch, bunch.copy()], opmd_diag=False)
 
 
 
diff --git a/tests/test_fluid_model.py b/tests/test_fluid_model.py
index d06fd3a..429115c 100644
--- a/tests/test_fluid_model.py
+++ b/tests/test_fluid_model.py
@@ -42,7 +42,7 @@ def test_fluid_model(plot=False):
 
     # Check final parameters.
     ene_sp = params_evolution['rel_ene_spread'][-1]
-    assert approx(ene_sp, rel=1e-10) == 0.024179998095119972
+    assert approx(ene_sp, rel=1e-10) == 0.024157374564016194
 
     # Quick plot of results.
     if plot:
diff --git a/tests/test_laser_envelope_model.py b/tests/test_laser_envelope_model.py
index addeefc..45648d5 100644
--- a/tests/test_laser_envelope_model.py
+++ b/tests/test_laser_envelope_model.py
@@ -3,6 +3,7 @@
 import scipy.constants as ct
 import matplotlib.pyplot as plt
 from wake_t.physics_models.laser.laser_pulse import GaussianPulse
+from wake_t.physics_models.laser.utils import unwrap
 
 
 def test_gaussian_laser_in_vacuum(plot=False):
@@ -190,6 +191,15 @@ def test_gaussian_laser_in_vacuum_with_subgrid(plot=False):
         plt.show()
 
 
+def test_unwrap():
+    """Test that the custom function for phase unwrapping implemented for
+    numba agrees with the numpy version.
+    """
+    phase = np.linspace(0, np.pi, num=5)
+    phase[3:] += np.pi
+    np.testing.assert_allclose(np.unwrap(phase), unwrap(phase))
+
+
 def calculate_spot_size(a_env, dr):
     # Project envelope to r
     a_proj = np.sum(np.abs(a_env), axis=0)
@@ -230,3 +240,4 @@ def calculate_a0(a_env):
 if __name__ == "__main__":
     test_gaussian_laser_in_vacuum(plot=True)
     test_gaussian_laser_in_vacuum_with_subgrid(plot=True)
+    test_unwrap()
diff --git a/tests/test_multibunch.py b/tests/test_multibunch.py
index 781837e..adb9bf6 100644
--- a/tests/test_multibunch.py
+++ b/tests/test_multibunch.py
@@ -55,8 +55,8 @@ def test_multibunch_plasma_simulation(plot=False):
     # Assert final parameters are correct.
     final_energy_driver = driver_params['avg_ene'][-1]
     final_energy_witness = witness_params['avg_ene'][-1]
-    assert approx(final_energy_driver, rel=1e-10) == 1700.3927190416732
-    assert approx(final_energy_witness, rel=1e-10) == 636.330857261392
+    assert approx(final_energy_driver, rel=1e-10) == 1700.3843657635728
+    assert approx(final_energy_witness, rel=1e-10) == 636.3260426124102
 
     if plot:
         z = driver_params['prop_dist'] * 1e2
diff --git a/tests/test_ramps.py b/tests/test_ramps.py
index 8e556ad..121723c 100644
--- a/tests/test_ramps.py
+++ b/tests/test_ramps.py
@@ -33,7 +33,7 @@ def test_downramp():
     downramp.track(bunch)
     bunch_params = analyze_bunch(bunch)
     beta_x = bunch_params['beta_x']
-    assert beta_x == 0.009750309290619276
+    assert beta_x == 0.009750308724018872
 
 
 def test_upramp():
@@ -64,7 +64,7 @@ def test_upramp():
     downramp.track(bunch)
     bunch_params = analyze_bunch(bunch)
     beta_x = bunch_params['beta_x']
-    assert beta_x == 0.0007631600676104024
+    assert beta_x == 0.000763155045965493
 
 
 if __name__ == '__main__':
diff --git a/tests/test_simple_blowout.py b/tests/test_simple_blowout.py
index 78d341d..c51d20e 100644
--- a/tests/test_simple_blowout.py
+++ b/tests/test_simple_blowout.py
@@ -45,7 +45,7 @@ def test_simple_blowout_wakefield(make_plots=False):
 
     bunch_params = analyze_bunch(bunch)
     rel_ene_sp = bunch_params['rel_ene_spread']
-    assert approx(rel_ene_sp, rel=1e-10) == 0.3637648484576557
+    assert approx(rel_ene_sp, rel=1e-10) == 0.3637651769109033
 
     if make_plots:
         # Analyze bunch evolution.
diff --git a/tutorials/00_basic_tutorial.py b/tutorials/00_basic_tutorial.py
index 5455887..5ab9ef4 100644
--- a/tutorials/00_basic_tutorial.py
+++ b/tutorials/00_basic_tutorial.py
@@ -63,7 +63,6 @@
 # average energy of :math:`100 \ \mathrm{MeV}` with a :math:`1 \ \%` spread
 # and a total charge of :math:`30 \ \mathrm{pC}`.
 
-from copy import deepcopy
 from wake_t.utilities.bunch_generation import get_gaussian_bunch_from_size
 
 # Beam parameters.
@@ -82,7 +81,7 @@
     q_bunch, n_part, name='elec_bunch')
 
 # Store bunch copy (will be needed later).
-bunch_bkp = deepcopy(bunch)
+bunch_bkp = bunch.copy()
 
 # Show phase space.
 bunch.show()
@@ -148,7 +147,7 @@
 # :math:`0.4 \ cm`, :math:`0.6 \ cm`, :math:`0.8 \ cm` and :math:`1.0 \ cm`.
 
 # Get again the original distribution.
-bunch = deepcopy(bunch_bkp)
+bunch = bunch_bkp.copy()
 
 # Create a 1 cm drift with 5 outputs (one every 0.2 cm).
 drift = Drift(length=1e-2, n_out=5)
@@ -195,7 +194,7 @@
 
 
 # Get again the original distribution.
-bunch = deepcopy(bunch_bkp)
+bunch = bunch_bkp.copy()
 
 # Create a 1 cm drift with 5 outputs (one every 0.2 cm).
 drift = Drift(length=1e-2, n_out=5)
diff --git a/tutorials/01_single_plasma_simulation.py b/tutorials/01_single_plasma_simulation.py
index 25f256d..2588f31 100644
--- a/tutorials/01_single_plasma_simulation.py
+++ b/tutorials/01_single_plasma_simulation.py
@@ -27,7 +27,6 @@
 # As a first step, let's generate a gaussian electron beam and keep a copy
 # of it for later use:
 
-from copy import deepcopy
 from wake_t.utilities.bunch_generation import get_gaussian_bunch_from_size
 
 # Beam parameters.
@@ -46,7 +45,7 @@
     q_bunch, n_part, name='elec_bunch')
 
 # Store bunch copy (will be needed later).
-bunch_bkp = deepcopy(bunch)
+bunch_bkp = bunch.copy()
 
 # Show phase space.
 bunch.show()
@@ -105,7 +104,7 @@
 from wake_t import GaussianPulse
 
 # Get again the original distribution.
-bunch = deepcopy(bunch_bkp)
+bunch = bunch_bkp.copy()
 
 # Laser parameters.
 laser_xi_c = 60e-6  # m (laser centroid in simulation box)
@@ -211,7 +210,7 @@ def density_profile(z):
 import scipy.constants as ct
 
 # Get again the original distribution.
-bunch = deepcopy(bunch_bkp)
+bunch = bunch_bkp.copy()
 
 # Calculate transverse parabolic profile.
 r_e = ct.e**2 / (4. * np.pi * ct.epsilon_0 * ct.m_e * ct.c**2)  # elec. radius
diff --git a/wake_t/beamline_elements/active_plasma_lens.py b/wake_t/beamline_elements/active_plasma_lens.py
index fa48f18..3647276 100644
--- a/wake_t/beamline_elements/active_plasma_lens.py
+++ b/wake_t/beamline_elements/active_plasma_lens.py
@@ -51,6 +51,20 @@ class ActivePlasmaLens(PlasmaStage):
         A list of values can also be provided. In this case, the list
         should have the same order as the list of bunches given to the
         ``track`` method.
+    push_bunches_before_diags : bool, optional
+        Whether to push the bunches before saving them to the diagnostics.
+        Since the time step of the diagnostics can be different from that
+        of the bunches, it could happen that the bunches appear in the
+        diagnostics as they were at the last push, but not at the actual
+        time of the diagnostics. Setting this parameter to ``True``
+        (default) ensures that an additional push is given to all bunches
+        to evolve them to the diagnostics time before saving.
+        This additional push will always have a time step smaller than
+        the the time step of the bunch, so it has no detrimental impact
+        on the accuracy of the simulation. However, it could make
+        convergence studies more difficult to interpret,
+        since the number of pushes will depend on `n_diags`. Therefore,
+        it is exposed as an option so that it can be disabled if needed.
     n_out : int
         Number of times along the lens in which the particle distribution
         should be returned (A list with all output bunches is returned
@@ -78,6 +92,7 @@ def __init__(
         wakefield_model: Optional[str] = 'quasistatic_2d',
         bunch_pusher: Optional[Literal['boris', 'rk4']] = 'boris',
         dt_bunch: Optional[DtBunchType] = 'auto',
+        push_bunches_before_diags: Optional[bool] = True,
         n_out: Optional[int] = 1,
         name: Optional[str] = 'Active plasma lens',
         **model_params
@@ -100,6 +115,7 @@ def __init__(
             wakefield_model=wakefield_model,
             bunch_pusher=bunch_pusher,
             dt_bunch=dt_bunch,
+            push_bunches_before_diags=push_bunches_before_diags,
             n_out=n_out,
             name=name,
             external_fields=[self.apl_field],
diff --git a/wake_t/beamline_elements/beamline.py b/wake_t/beamline_elements/beamline.py
index b7dc4b3..eb37b2f 100644
--- a/wake_t/beamline_elements/beamline.py
+++ b/wake_t/beamline_elements/beamline.py
@@ -20,7 +20,8 @@ def track(
         self,
         bunches: Optional[Union[ParticleBunch, List[ParticleBunch]]] = [],
         opmd_diag: Optional[bool] = False,
-        diag_dir: Optional[str] = None
+        diag_dir: Optional[str] = None,
+        show_progress_bar: Optional[bool] = True,
     ) -> Union[List[ParticleBunch], List[List[ParticleBunch]]]:
         """
         Track bunch through beamline.
@@ -40,6 +41,9 @@ def track(
             Directory into which the openPMD output will be written. By default
             this is a 'diags' folder in the current directory. Only needed if
             `opmd_diag=True`.
+        show_progress_bar : bool, optional
+            Whether to show a progress bar of the tracking through each
+            element. By default ``True``.
 
         Returns
         -------
@@ -50,5 +54,11 @@ def track(
         if type(opmd_diag) is not OpenPMDDiagnostics and opmd_diag:
             opmd_diag = OpenPMDDiagnostics(write_dir=diag_dir)
         for element in self.elements:
-            bunch_list.extend(element.track(bunches, opmd_diag=opmd_diag))
+            bunch_list.extend(
+                element.track(
+                    bunches,
+                    opmd_diag=opmd_diag,
+                    show_progress_bar=show_progress_bar,
+                )
+            )
         return bunch_list
diff --git a/wake_t/beamline_elements/field_element.py b/wake_t/beamline_elements/field_element.py
index 15fe6c1..31159d9 100644
--- a/wake_t/beamline_elements/field_element.py
+++ b/wake_t/beamline_elements/field_element.py
@@ -39,6 +39,20 @@ class FieldElement():
         Function used to determine the adaptive time step for bunches in
         which the time step is set to ``'auto'``. The function should take
         solely a ``ParticleBunch`` as argument.
+    push_bunches_before_diags : bool, optional
+        Whether to push the bunches before saving them to the diagnostics.
+        Since the time step of the diagnostics can be different from that
+        of the bunches, it could happen that the bunches appear in the
+        diagnostics as they were at the last push, but not at the actual
+        time of the diagnostics. Setting this parameter to ``True``
+        (default) ensures that an additional push is given to all bunches
+        to evolve them to the diagnostics time before saving.
+        This additional push will always have a time step smaller than
+        the the time step of the bunch, so it has no detrimental impact
+        on the accuracy of the simulation. However, it could make
+        convergence studies more difficult to interpret,
+        since the number of pushes will depend on `n_diags`. Therefore,
+        it is exposed as an option so that it can be disabled if needed.
 
     """
 
@@ -50,7 +64,8 @@ def __init__(
         n_out: Optional[int] = 1,
         name: Optional[str] = 'field element',
         fields: Optional[List[Field]] = [],
-        auto_dt_bunch: Optional[str] = None
+        auto_dt_bunch: Optional[str] = None,
+        push_bunches_before_diags: Optional[bool] = True,
     ) -> None:
         self.length = length
         self.bunch_pusher = bunch_pusher
@@ -59,12 +74,14 @@ def __init__(
         self.name = name
         self.fields = fields
         self.auto_dt_bunch = auto_dt_bunch
+        self.push_bunches_before_diags = push_bunches_before_diags
 
     def track(
         self,
         bunches: Optional[Union[ParticleBunch, List[ParticleBunch]]] = [],
         opmd_diag: Optional[Union[bool, OpenPMDDiagnostics]] = False,
-        diag_dir: Optional[str] = None
+        diag_dir: Optional[str] = None,
+        show_progress_bar: Optional[bool] = True,
     ) -> Union[List[ParticleBunch], List[List[ParticleBunch]]]:
         """
         Track bunch through element.
@@ -83,6 +100,9 @@ def track(
             Directory into which the openPMD output will be written. By default
             this is a 'diags' folder in the current directory. Only needed if
             `opmd_diag=True`.
+        show_progress_bar : bool, optional
+            Whether to show a progress bar of the tracking. By default
+            ``True``.
 
         Returns
         -------
@@ -119,6 +139,8 @@ def track(
             opmd_diags=opmd_diag,
             bunch_pusher=self.bunch_pusher,
             auto_dt_bunch_f=self.auto_dt_bunch,
+            push_bunches_before_diags=self.push_bunches_before_diags,
+            show_progress_bar=show_progress_bar,
             section_name=self.name
         )
 
diff --git a/wake_t/beamline_elements/plasma_ramp.py b/wake_t/beamline_elements/plasma_ramp.py
index ed3376a..3590165 100644
--- a/wake_t/beamline_elements/plasma_ramp.py
+++ b/wake_t/beamline_elements/plasma_ramp.py
@@ -96,6 +96,20 @@ class PlasmaRamp(PlasmaStage):
         A list of values can also be provided. In this case, the list
         should have the same order as the list of bunches given to the
         ``track`` method.
+    push_bunches_before_diags : bool, optional
+        Whether to push the bunches before saving them to the diagnostics.
+        Since the time step of the diagnostics can be different from that
+        of the bunches, it could happen that the bunches appear in the
+        diagnostics as they were at the last push, but not at the actual
+        time of the diagnostics. Setting this parameter to ``True``
+        (default) ensures that an additional push is given to all bunches
+        to evolve them to the diagnostics time before saving.
+        This additional push will always have a time step smaller than
+        the the time step of the bunch, so it has no detrimental impact
+        on the accuracy of the simulation. However, it could make
+        convergence studies more difficult to interpret,
+        since the number of pushes will depend on `n_diags`. Therefore,
+        it is exposed as an option so that it can be disabled if needed.
     n_out : int
         Number of times along the stage in which the particle distribution
         should be returned (A list with all output bunches is returned
@@ -126,6 +140,7 @@ def __init__(
         position_down: Optional[float] = None,
         bunch_pusher: Optional[Literal['boris', 'rk4']] = 'boris',
         dt_bunch: Optional[DtBunchType] = 'auto',
+        push_bunches_before_diags: Optional[bool] = True,
         n_out: Optional[int] = 1,
         name: Optional[str] = 'Plasma ramp',
         **model_params
@@ -151,6 +166,7 @@ def __init__(
             wakefield_model=wakefield_model,
             bunch_pusher=bunch_pusher,
             dt_bunch=dt_bunch,
+            push_bunches_before_diags=push_bunches_before_diags,
             n_out=n_out,
             name=name,
             **model_params
diff --git a/wake_t/beamline_elements/plasma_stage.py b/wake_t/beamline_elements/plasma_stage.py
index 69d6420..ec9a83f 100644
--- a/wake_t/beamline_elements/plasma_stage.py
+++ b/wake_t/beamline_elements/plasma_stage.py
@@ -48,6 +48,20 @@ class PlasmaStage(FieldElement):
         stage. A list of values can also be provided. In this case, the list
         should have the same order as the list of bunches given to the
         ``track`` method.
+    push_bunches_before_diags : bool, optional
+        Whether to push the bunches before saving them to the diagnostics.
+        Since the time step of the diagnostics can be different from that
+        of the bunches, it could happen that the bunches appear in the
+        diagnostics as they were at the last push, but not at the actual
+        time of the diagnostics. Setting this parameter to ``True``
+        (default) ensures that an additional push is given to all bunches
+        to evolve them to the diagnostics time before saving.
+        This additional push will always have a time step smaller than
+        the the time step of the bunch, so it has no detrimental impact
+        on the accuracy of the simulation. However, it could make
+        convergence studies more difficult to interpret,
+        since the number of pushes will depend on `n_diags`. Therefore,
+        it is exposed as an option so that it can be disabled if needed.
     n_out : int, optional
         Number of times along the stage in which the particle distribution
         should be returned (A list with all output bunches is returned
@@ -77,6 +91,7 @@ def __init__(
         wakefield_model: Optional[str] = 'simple_blowout',
         bunch_pusher: Optional[Literal['boris', 'rk4']] = 'boris',
         dt_bunch: Optional[DtBunchType] = 'auto',
+        push_bunches_before_diags: Optional[bool] = True,
         n_out: Optional[int] = 1,
         name: Optional[str] = 'Plasma stage',
         external_fields: Optional[List[Field]] = [],
@@ -96,7 +111,8 @@ def __init__(
             n_out=n_out,
             name=name,
             fields=fields,
-            auto_dt_bunch=self._get_optimized_dt
+            auto_dt_bunch=self._get_optimized_dt,
+            push_bunches_before_diags=push_bunches_before_diags,
         )
 
     def _get_density_profile(self, density):
diff --git a/wake_t/beamline_elements/tm_elements.py b/wake_t/beamline_elements/tm_elements.py
index 843d5bb..bf1647d 100644
--- a/wake_t/beamline_elements/tm_elements.py
+++ b/wake_t/beamline_elements/tm_elements.py
@@ -1,7 +1,6 @@
 """ Contains the classes of all elements tracked using transfer matrices. """
 from typing import Optional, Union, List
 import time
-from copy import deepcopy
 
 import numpy as np
 import scipy.constants as ct
@@ -77,7 +76,8 @@ def track(
         backtrack: Optional[bool] = False,
         out_initial: Optional[bool] = False,
         opmd_diag: Optional[Union[bool, OpenPMDDiagnostics]] = False,
-        diag_dir: Optional[str] = None
+        diag_dir: Optional[str] = None,
+        show_progress_bar: Optional[bool] = True,
     ) -> List[ParticleBunch]:
         """
         Track bunch through element.
@@ -100,6 +100,9 @@ def track(
             Directory into which the openPMD output will be written. By default
             this is a 'diags' folder in the current directory. Only needed if
             `opmd_diag=True`.
+        show_progress_bar : bool, optional
+            Whether to show a progress bar of the tracking. By default
+            ``True``.
 
         Returns
         -------
@@ -126,26 +129,28 @@ def track(
             opmd_diag = None
 
         # Print output header
-        print('')
-        print(self.element_name.capitalize())
-        print('-'*len(self.element_name))
-        self._print_element_properties()
-        csr_string = 'on' if self.csr_on else 'off'
-        print('CSR {}.'.format(csr_string))
-        print('')
-        n_steps = len(track_steps)
-        st_0 = 'Tracking in {} step(s)... '.format(n_steps)
+        if show_progress_bar:
+            print('')
+            print(self.element_name.capitalize())
+            print('-'*len(self.element_name))
+            self._print_element_properties()
+            csr_string = 'on' if self.csr_on else 'off'
+            print('CSR {}.'.format(csr_string))
+            print('')
+            n_steps = len(track_steps)
+            st_0 = 'Tracking in {} step(s)... '.format(n_steps)
 
         # Start tracking
         start_time = time.time()
         output_bunch_list = list()
         if out_initial:
-            output_bunch_list.append(deepcopy(bunch))
+            output_bunch_list.append(bunch.copy())
             if opmd_diag is not None:
                 opmd_diag.write_diagnostics(
                     0., l_step/ct.c, [output_bunch_list[-1]])
         for i in track_steps:
-            print_progress_bar(st_0, i+1, n_steps)
+            if show_progress_bar:
+                print_progress_bar(st_0, i+1, n_steps)
             l_curr = (i+1) * l_step * (1-2*backtrack)
             # Track with transfer matrix
             bunch_mat = track_with_transfer_map(
@@ -174,9 +179,10 @@ def track(
             opmd_diag.increase_z_pos(self.length)
 
         # Finalize
-        tracking_time = time.time() - start_time
-        print('Done ({} s).'.format(tracking_time))
-        print('-'*80)
+        if show_progress_bar:
+            tracking_time = time.time() - start_time
+            print('Done ({} s).'.format(tracking_time))
+            print('-'*80)
         return output_bunch_list
 
     def _get_beam_matrix_for_tracking(self, bunch):
@@ -219,7 +225,7 @@ def _update_input_bunch(self, bunch, bunch_mat, output_bunch_list):
             last_bunch = self._create_new_bunch(bunch, new_bunch_mat,
                                                 self.length)
         else:
-            last_bunch = deepcopy(output_bunch_list[-1])
+            last_bunch = output_bunch_list[-1].copy()
         bunch.set_phase_space(last_bunch.x, last_bunch.y, last_bunch.xi,
                               last_bunch.px, last_bunch.py, last_bunch.pz)
         bunch.prop_distance = last_bunch.prop_distance
diff --git a/wake_t/particles/particle_bunch.py b/wake_t/particles/particle_bunch.py
index c212e10..09e3018 100644
--- a/wake_t/particles/particle_bunch.py
+++ b/wake_t/particles/particle_bunch.py
@@ -2,6 +2,8 @@
 This module contains the class defining a particle bunch.
 """
 # TODO: clean methods to set and get bunch matrix
+from __future__ import annotations
+from copy import deepcopy
 from typing import Optional
 
 import numpy as np
@@ -295,6 +297,29 @@ def evolve(self, fields, t, dt, pusher='rk4'):
             )
         self.prop_distance += dt * ct.c
 
+    def copy(self) -> ParticleBunch:
+        """Return a copy of the bunch.
+
+        To improve performance, this copy won't contain copies of auxiliary
+        arrays, only of the particle coordinates and properties.
+        """
+        bunch_copy = ParticleBunch(
+            w=deepcopy(self.w),
+            x=deepcopy(self.x),
+            y=deepcopy(self.y),
+            xi=deepcopy(self.xi),
+            px=deepcopy(self.px),
+            py=deepcopy(self.py),
+            pz=deepcopy(self.pz),
+            prop_distance=deepcopy(self.prop_distance),
+            name=deepcopy(self.name),
+            q_species=deepcopy(self.q_species),
+            m_species=deepcopy(self.m_species)
+        )
+        bunch_copy.x_ref = self.x_ref
+        bunch_copy.theta_ref = self.theta_ref
+        return bunch_copy
+
     def get_field_arrays(self):
         """Get the arrays where the gathered fields will be stored."""
         if not self.__field_arrays_allocated:
diff --git a/wake_t/tracking/progress_bar.py b/wake_t/tracking/progress_bar.py
index 8ea08df..d050a22 100644
--- a/wake_t/tracking/progress_bar.py
+++ b/wake_t/tracking/progress_bar.py
@@ -10,7 +10,7 @@
 warnings.filterwarnings('ignore', '.*clamping.*', )
 
 
-def get_progress_bar(description, total_length):
+def get_progress_bar(description, total_length, disable):
     """Get progress bar for the tracker.
 
     Parameters
@@ -19,6 +19,8 @@ def get_progress_bar(description, total_length):
         Description to be appended to start of the progress bar.
     total_length : float
         Total length in metres of the stage to be tracked.
+    disable : bool
+        Whether to disable (not show) the progress bar.
 
     Returns
     -------
@@ -31,6 +33,7 @@ def get_progress_bar(description, total_length):
         total=total_length,
         unit='m',
         bar_format=l_bar + "{bar}" + r_bar,
-        file=sys.stdout
+        file=sys.stdout,
+        disable=disable
     )
     return progress_bar
diff --git a/wake_t/tracking/tracker.py b/wake_t/tracking/tracker.py
index b9dcd69..0a3600e 100644
--- a/wake_t/tracking/tracker.py
+++ b/wake_t/tracking/tracker.py
@@ -1,6 +1,5 @@
 """ This module contains the Tracker class. """
 from typing import Optional, Callable, List, Literal
-from copy import deepcopy
 
 import numpy as np
 import scipy.constants as ct
@@ -62,6 +61,22 @@ class Tracker():
     bunch_pusher : str, optional
         The particle pusher used to evolve the bunches. Possible values
         are `'boris'` or `'rk4'`.
+    push_bunches_before_diags : bool, optional
+        Whether to push the bunches before saving them to the diagnostics.
+        Since the time step of the diagnostics can be different from that
+        of the bunches, it could happen that the bunches appear in the
+        diagnostics as they were at the last push, but not at the actual
+        time of the diagnostics. Setting this parameter to ``True``
+        (default) ensures that an additional push is given to all bunches
+        to evolve them to the diagnostics time before saving.
+        This additional push will always have a time step smaller than
+        the the time step of the bunch, so it has no detrimental impact
+        on the accuracy of the simulation. However, it could make
+        convergence studies more difficult to interpret,
+        since the number of pushes will depend on `n_diags`. Therefore,
+        it is exposed as an option so that it can be disabled if needed.
+    show_progress_bar : bool, optional
+        Whether to show a progress bar of the tracking. By default ``True``.
     section_name : str, optional
         Name of the section to be tracked. This will be appended to the
         beginning of the progress bar.
@@ -78,6 +93,8 @@ def __init__(
         opmd_diags: Optional[OpenPMDDiagnostics] = None,
         auto_dt_bunch_f: Optional[Callable[[ParticleBunch], float]] = None,
         bunch_pusher: Optional[Literal['boris', 'rk4']] = 'boris',
+        push_bunches_before_diags: Optional[bool] = True,
+        show_progress_bar: Optional[bool] = True,
         section_name: Optional[str] = 'Simulation'
     ) -> None:
         self.t_final = t_final
@@ -88,6 +105,8 @@ def __init__(
         self.n_diags = n_diags
         self.auto_dt_bunch_f = auto_dt_bunch_f
         self.bunch_pusher = bunch_pusher
+        self.push_bunches_before_diags = push_bunches_before_diags
+        self.show_progress_bar = show_progress_bar
         self.section_name = section_name
 
         # Get all numerical fields and their time steps.
@@ -134,7 +153,11 @@ def do_tracking(self) -> List[List[ParticleBunch]]:
         set_num_threads(num_threads)
 
         # Initialize progress bar.
-        progress_bar = get_progress_bar(self.section_name, self.t_final*ct.c)
+        progress_bar = get_progress_bar(
+            description=self.section_name,
+            total_length=self.t_final*ct.c,
+            disable=not self.show_progress_bar,
+        )
 
         # Calculate fields at t=0.
         for field in self.num_fields:
@@ -187,20 +210,14 @@ def do_tracking(self) -> List[List[ParticleBunch]]:
 
             # If next object is a ParticleBunch, update it.
             if isinstance(obj_next, ParticleBunch):
-                obj_next.evolve(
-                    self.fields, t_current, dt_next, self.bunch_pusher)
-                # Update the time step if set to `'auto'`.
-                if obj_next in self.auto_dt_bunches:
-                    dt_objects[i_next] = self.auto_dt_bunch_f(obj_next)
-                # Determine if this was the last push.
-                final_push = np.float32(t_next) == np.float32(self.t_final)
-                # Determine if next push brings the bunch beyond `t_final`.
-                next_push_beyond_final_time = (
-                    t_next + dt_objects[i_next] > self.t_final)
-                # Make sure the last push of the bunch advances it to exactly
-                # `t_final`.
-                if not final_push and next_push_beyond_final_time:
-                    dt_objects[i_next] = self.t_final - t_next
+                self.evolve_bunch(
+                    bunch=obj_next,
+                    t_current=t_current,
+                    t_next=t_next,
+                    dt_next=dt_next,
+                    i_next=i_next,
+                    dt_objects=dt_objects,
+                )
 
             # If next object is a NumericalField, update it.
             elif isinstance(obj_next, NumericalField):
@@ -208,6 +225,20 @@ def do_tracking(self) -> List[List[ParticleBunch]]:
 
             # If next object are the diagnostics, generate them.
             elif obj_next == 'diags':
+                # Evolve all bunches to the diagnostics time.
+                if self.push_bunches_before_diags:
+                    for i, obj in enumerate(self.objects_to_track):
+                        if isinstance(obj, ParticleBunch):
+                            dt_bunch = t_next - t_objects[i]
+                            self.evolve_bunch(
+                                bunch=obj,
+                                t_current=t_objects[i],
+                                t_next=t_next,
+                                dt_next=dt_bunch,
+                                i_next=i,
+                                dt_objects=dt_objects,
+                            )
+                            t_objects[i] += dt_bunch
                 self.generate_diagnostics()
 
             # Advance current time of the update object.
@@ -230,25 +261,51 @@ def do_tracking(self) -> List[List[ParticleBunch]]:
 
         return self.bunch_list
 
+    def evolve_bunch(
+        self,
+        bunch: ParticleBunch,
+        t_current: float,
+        t_next: float,
+        dt_next: float,
+        i_next: int,
+        dt_objects: List,
+    ):
+        """Evolve particle bunch to next time step.
+
+        Parameters
+        ----------
+        bunch : ParticleBunch
+            The particle bunch to evolve.
+        t_current : float
+            The current time of the simulation.
+        t_next : float
+            The time of the next step.
+        dt_next : float
+            The time step by which to advance the bunch.
+        i_next : int
+            The index of the bunch in the list of objects to track.
+        dt_objects : List
+            The time steps of all objects to track.
+        """
+        bunch.evolve(self.fields, t_current, dt_next, self.bunch_pusher)
+        # Update the time step if set to `'auto'`.
+        if bunch in self.auto_dt_bunches:
+            dt_objects[i_next] = self.auto_dt_bunch_f(bunch)
+        # Determine if this was the last push.
+        final_push = np.float32(t_next) == np.float32(self.t_final)
+        # Determine if next push brings the bunch beyond `t_final`.
+        next_push_beyond_final_time = (
+            t_next + dt_objects[i_next] > self.t_final)
+        # Make sure the last push of the bunch advances it to exactly
+        # `t_final`.
+        if not final_push and next_push_beyond_final_time:
+            dt_objects[i_next] = self.t_final - t_next
+
     def generate_diagnostics(self) -> None:
         """Generate tracking diagnostics."""
         # Make copy of current bunches and store in output list.
         for i, bunch in enumerate(self.bunches):
-            self.bunch_list[i].append(
-                ParticleBunch(
-                    deepcopy(bunch.w),
-                    deepcopy(bunch.x),
-                    deepcopy(bunch.y),
-                    deepcopy(bunch.xi),
-                    deepcopy(bunch.px),
-                    deepcopy(bunch.py),
-                    deepcopy(bunch.pz),
-                    prop_distance=deepcopy(bunch.prop_distance),
-                    name=deepcopy(bunch.name),
-                    q_species=deepcopy(bunch.q_species),
-                    m_species=deepcopy(bunch.m_species)
-                )
-            )
+            self.bunch_list[i].append(bunch.copy())
 
         # If needed, write also the openPMD diagnostics.
         if self.opmd_diags is not None: