diff --git a/cirq-core/cirq/circuits/circuit_operation_test.py b/cirq-core/cirq/circuits/circuit_operation_test.py
index c0596ab261c..e708202c46c 100644
--- a/cirq-core/cirq/circuits/circuit_operation_test.py
+++ b/cirq-core/cirq/circuits/circuit_operation_test.py
@@ -267,7 +267,7 @@ def test_recursive_params():
 
     # First example should behave like an X when simulated
     result = cirq.Simulator().simulate(cirq.Circuit(circuitop), param_resolver=outer_params)
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
 
 
 @pytest.mark.parametrize('add_measurements', [True, False])
@@ -343,9 +343,9 @@ def test_repeat_zero_times(add_measurements, use_repetition_ids, initial_reps):
         subcircuit.freeze(), repetitions=initial_reps, use_repetition_ids=use_repetition_ids
     )
     result = cirq.Simulator().simulate(cirq.Circuit(op))
-    assert np.allclose(result.state_vector(), [0, 1] if initial_reps % 2 else [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [0, 1] if initial_reps % 2 else [1, 0])
     result = cirq.Simulator().simulate(cirq.Circuit(op**0))
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
 
 
 def test_no_repetition_ids():
@@ -375,13 +375,13 @@ def test_parameterized_repeat():
     assert cirq.parameter_names(op) == {'a'}
     assert not cirq.has_unitary(op)
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 0})
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1})
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 2})
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': -1})
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
     with pytest.raises(TypeError, match='Only integer or sympy repetitions are allowed'):
         cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1.5})
     with pytest.raises(ValueError, match='Circuit contains ops whose symbols were not specified'):
@@ -390,13 +390,13 @@ def test_parameterized_repeat():
     assert cirq.parameter_names(op) == {'a'}
     assert not cirq.has_unitary(op)
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 0})
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1})
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 2})
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': -1})
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
     with pytest.raises(TypeError, match='Only integer or sympy repetitions are allowed'):
         cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1.5})
     with pytest.raises(ValueError, match='Circuit contains ops whose symbols were not specified'):
@@ -405,11 +405,11 @@ def test_parameterized_repeat():
     assert cirq.parameter_names(op) == {'a', 'b'}
     assert not cirq.has_unitary(op)
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1, 'b': 1})
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 2, 'b': 1})
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1, 'b': 2})
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
     with pytest.raises(TypeError, match='Only integer or sympy repetitions are allowed'):
         cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1.5, 'b': 1})
     with pytest.raises(ValueError, match='Circuit contains ops whose symbols were not specified'):
@@ -418,11 +418,11 @@ def test_parameterized_repeat():
     assert cirq.parameter_names(op) == {'a', 'b'}
     assert not cirq.has_unitary(op)
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1, 'b': 1})
-    assert np.allclose(result.state_vector(), [1, 0])
+    assert np.allclose(result.state_vector(copy=False), [1, 0])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1.5, 'b': 1})
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
     result = cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1, 'b': 1.5})
-    assert np.allclose(result.state_vector(), [0, 1])
+    assert np.allclose(result.state_vector(copy=False), [0, 1])
     with pytest.raises(TypeError, match='Only integer or sympy repetitions are allowed'):
         cirq.Simulator().simulate(cirq.Circuit(op), param_resolver={'a': 1.5, 'b': 1.5})
     with pytest.raises(ValueError, match='Circuit contains ops whose symbols were not specified'):
diff --git a/cirq-core/cirq/contrib/quantum_volume/quantum_volume.py b/cirq-core/cirq/contrib/quantum_volume/quantum_volume.py
index 042a895d4d4..9dd947727e5 100644
--- a/cirq-core/cirq/contrib/quantum_volume/quantum_volume.py
+++ b/cirq-core/cirq/contrib/quantum_volume/quantum_volume.py
@@ -81,12 +81,14 @@ def compute_heavy_set(circuit: cirq.Circuit) -> List[int]:
     # output is defined in terms of probabilities, where our wave function is in
     # terms of amplitudes. We convert it by using the Born rule: squaring each
     # amplitude and taking their absolute value
-    median = np.median(np.abs(results.state_vector() ** 2))
+    median = np.median(np.abs(results.state_vector(copy=False) ** 2))
 
     # The output wave function is a vector from the result value (big-endian) to
     # the probability of that bit-string. Return all of the bit-string
     # values that have a probability greater than the median.
-    return [idx for idx, amp in enumerate(results.state_vector()) if np.abs(amp**2) > median]
+    return [
+        idx for idx, amp in enumerate(results.state_vector(copy=False)) if np.abs(amp**2) > median
+    ]
 
 
 @dataclass
diff --git a/cirq-core/cirq/experiments/grid_parallel_two_qubit_xeb.py b/cirq-core/cirq/experiments/grid_parallel_two_qubit_xeb.py
index d1bb1283383..8c1b2d2674f 100644
--- a/cirq-core/cirq/experiments/grid_parallel_two_qubit_xeb.py
+++ b/cirq-core/cirq/experiments/grid_parallel_two_qubit_xeb.py
@@ -504,7 +504,8 @@ def _get_xeb_result(
             while moment_index < 2 * depth:
                 step_result = next(step_results)
                 moment_index += 1
-            amplitudes = step_result.state_vector()
+            # copy=False is safe because state_vector_to_probabilities will copy anyways
+            amplitudes = step_result.state_vector(copy=False)
             probabilities = value.state_vector_to_probabilities(amplitudes)
             _, counts = np.unique(measurements, return_counts=True)
             empirical_probs = counts / len(measurements)
diff --git a/cirq-core/cirq/experiments/xeb_simulation.py b/cirq-core/cirq/experiments/xeb_simulation.py
index 72bbf252217..7fd19c530a7 100644
--- a/cirq-core/cirq/experiments/xeb_simulation.py
+++ b/cirq-core/cirq/experiments/xeb_simulation.py
@@ -65,7 +65,8 @@ def __call__(self, task: _Simulate2qXEBTask) -> List[Dict[str, Any]]:
             if cycle_depth not in cycle_depths:
                 continue
 
-            psi = step_result.state_vector()
+            # copy=False is safe because state_vector_to_probabilities will copy anyways
+            psi = step_result.state_vector(copy=False)
             pure_probs = value.state_vector_to_probabilities(psi)
 
             records += [
diff --git a/cirq-core/cirq/ops/boolean_hamiltonian_test.py b/cirq-core/cirq/ops/boolean_hamiltonian_test.py
index c4a3f781720..426f7a6b9eb 100644
--- a/cirq-core/cirq/ops/boolean_hamiltonian_test.py
+++ b/cirq-core/cirq/ops/boolean_hamiltonian_test.py
@@ -80,7 +80,11 @@ def test_circuit(boolean_str):
 
     circuit.append(hamiltonian_gate.on(*qubits))
 
-    phi = cirq.Simulator().simulate(circuit, qubit_order=qubits, initial_state=0).state_vector()
+    phi = (
+        cirq.Simulator()
+        .simulate(circuit, qubit_order=qubits, initial_state=0)
+        .state_vector(copy=False)
+    )
     actual = np.arctan2(phi.real, phi.imag) - math.pi / 2.0 > 0.0
 
     # Compare the two:
diff --git a/cirq-core/cirq/ops/common_gates_test.py b/cirq-core/cirq/ops/common_gates_test.py
index f268fdf6a27..f62ef507d17 100644
--- a/cirq-core/cirq/ops/common_gates_test.py
+++ b/cirq-core/cirq/ops/common_gates_test.py
@@ -1088,7 +1088,7 @@ def test_xpow_dim_3():
 
     sim = cirq.Simulator()
     circuit = cirq.Circuit([x(cirq.LineQid(0, 3)) ** 0.5] * 6)
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit)]
+    svs = [step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit)]
     # fmt: off
     expected = [
         [0.67, 0.67, 0.33],
@@ -1116,7 +1116,7 @@ def test_xpow_dim_4():
 
     sim = cirq.Simulator()
     circuit = cirq.Circuit([x(cirq.LineQid(0, 4)) ** 0.5] * 8)
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit)]
+    svs = [step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit)]
     # fmt: off
     expected = [
         [0.65, 0.65, 0.27, 0.27],
@@ -1147,11 +1147,15 @@ def test_zpow_dim_3():
 
     sim = cirq.Simulator()
     circuit = cirq.Circuit([z(cirq.LineQid(0, 3)) ** 0.5] * 6)
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit, initial_state=0)]
+    svs = [
+        step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit, initial_state=0)
+    ]
     expected = [[1, 0, 0]] * 6
     assert np.allclose((svs), expected)
 
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit, initial_state=1)]
+    svs = [
+        step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit, initial_state=1)
+    ]
     # fmt: off
     expected = [
         [0, L**0.5, 0],
@@ -1164,7 +1168,9 @@ def test_zpow_dim_3():
     # fmt: on
     assert np.allclose((svs), expected)
 
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit, initial_state=2)]
+    svs = [
+        step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit, initial_state=2)
+    ]
     # fmt: off
     expected = [
         [0, 0, L],
@@ -1192,11 +1198,15 @@ def test_zpow_dim_4():
 
     sim = cirq.Simulator()
     circuit = cirq.Circuit([z(cirq.LineQid(0, 4)) ** 0.5] * 8)
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit, initial_state=0)]
+    svs = [
+        step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit, initial_state=0)
+    ]
     expected = [[1, 0, 0, 0]] * 8
     assert np.allclose((svs), expected)
 
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit, initial_state=1)]
+    svs = [
+        step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit, initial_state=1)
+    ]
     # fmt: off
     expected = [
         [0, 1j**0.5, 0, 0],
@@ -1211,7 +1221,9 @@ def test_zpow_dim_4():
     # fmt: on
     assert np.allclose(svs, expected)
 
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit, initial_state=2)]
+    svs = [
+        step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit, initial_state=2)
+    ]
     # fmt: off
     expected = [
         [0, 0, 1j, 0],
@@ -1226,7 +1238,9 @@ def test_zpow_dim_4():
     # fmt: on
     assert np.allclose(svs, expected)
 
-    svs = [step.state_vector() for step in sim.simulate_moment_steps(circuit, initial_state=3)]
+    svs = [
+        step.state_vector(copy=True) for step in sim.simulate_moment_steps(circuit, initial_state=3)
+    ]
     # fmt: off
     expected = [
         [0, 0, 0, 1j**1.5],
diff --git a/cirq-core/cirq/sim/mux.py b/cirq-core/cirq/sim/mux.py
index f375dba580d..c0e357284de 100644
--- a/cirq-core/cirq/sim/mux.py
+++ b/cirq-core/cirq/sim/mux.py
@@ -168,7 +168,7 @@ def final_state_vector(
         param_resolver=param_resolver,
     )
 
-    return result.state_vector()
+    return result.state_vector(copy=False)
 
 
 def sample_sweep(
diff --git a/cirq-core/cirq/sim/simulator_test.py b/cirq-core/cirq/sim/simulator_test.py
index 68447b3ebe3..7640f4d228d 100644
--- a/cirq-core/cirq/sim/simulator_test.py
+++ b/cirq-core/cirq/sim/simulator_test.py
@@ -416,7 +416,9 @@ def _kraus_(self):
             cirq.Circuit(Reset11To00().on(*cirq.LineQubit.range(2))), initial_state=k
         )
         np.testing.assert_allclose(
-            out.state_vector(), cirq.one_hot(index=k % 3, shape=4, dtype=np.complex64), atol=1e-8
+            out.state_vector(copy=False),
+            cirq.one_hot(index=k % 3, shape=4, dtype=np.complex64),
+            atol=1e-8,
         )
 
 
diff --git a/cirq-core/cirq/sim/sparse_simulator.py b/cirq-core/cirq/sim/sparse_simulator.py
index 48eee8d33ae..d3c32d2022e 100644
--- a/cirq-core/cirq/sim/sparse_simulator.py
+++ b/cirq-core/cirq/sim/sparse_simulator.py
@@ -18,7 +18,7 @@
 
 import numpy as np
 
-from cirq import ops
+from cirq import _compat, ops
 from cirq._compat import deprecated_parameter
 from cirq.sim import simulator, state_vector, state_vector_simulator, state_vector_simulation_state
 
@@ -246,7 +246,7 @@ def __init__(
         self._dtype = dtype
         self._state_vector: Optional[np.ndarray] = None
 
-    def state_vector(self, copy: bool = True):
+    def state_vector(self, copy: Optional[bool] = None):
         """Return the state vector at this point in the computation.
 
         The state is returned in the computational basis with these basis
@@ -279,6 +279,12 @@ def state_vector(self, copy: bool = True):
                 parameters from the state vector and store then using False
                 can speed up simulation by eliminating a memory copy.
         """
+        if copy is None:
+            _compat._warn_or_error(
+                "Starting in v0.16, state_vector will not copy the state by default. "
+                "Explicitly set copy=True to copy the state."
+            )
+            copy = True
         if self._state_vector is None:
             self._state_vector = np.array([1])
             state = self._merged_sim_state
diff --git a/cirq-core/cirq/sim/sparse_simulator_test.py b/cirq-core/cirq/sim/sparse_simulator_test.py
index ff32871acd5..14966a4bdd2 100644
--- a/cirq-core/cirq/sim/sparse_simulator_test.py
+++ b/cirq-core/cirq/sim/sparse_simulator_test.py
@@ -550,9 +550,20 @@ def test_simulate_moment_steps(dtype: Type[np.number], split: bool):
     simulator = cirq.Simulator(dtype=dtype, split_untangled_states=split)
     for i, step in enumerate(simulator.simulate_moment_steps(circuit)):
         if i == 0:
-            np.testing.assert_almost_equal(step.state_vector(), np.array([0.5] * 4))
+            np.testing.assert_almost_equal(step.state_vector(copy=True), np.array([0.5] * 4))
         else:
-            np.testing.assert_almost_equal(step.state_vector(), np.array([1, 0, 0, 0]))
+            np.testing.assert_almost_equal(step.state_vector(copy=True), np.array([1, 0, 0, 0]))
+
+
+def test_simulate_moment_steps_implicit_copy_deprecated():
+    q0 = cirq.LineQubit(0)
+    simulator = cirq.Simulator()
+    steps = list(simulator.simulate_moment_steps(cirq.Circuit(cirq.X(q0))))
+
+    with cirq.testing.assert_deprecated(
+        "state_vector will not copy the state by default", deadline="v0.16"
+    ):
+        _ = steps[0].state_vector()
 
 
 @pytest.mark.parametrize('dtype', [np.complex64, np.complex128])
@@ -563,7 +574,7 @@ def test_simulate_moment_steps_empty_circuit(dtype: Type[np.number], split: bool
     step = None
     for step in simulator.simulate_moment_steps(circuit):
         pass
-    assert np.allclose(step.state_vector(), np.array([1]))
+    assert np.allclose(step.state_vector(copy=True), np.array([1]))
     assert not step.qubit_map
 
 
@@ -599,10 +610,10 @@ def test_simulate_moment_steps_intermediate_measurement(dtype: Type[np.number],
             result = int(step.measurements['q(0)'][0])
             expected = np.zeros(2)
             expected[result] = 1
-            np.testing.assert_almost_equal(step.state_vector(), expected)
+            np.testing.assert_almost_equal(step.state_vector(copy=True), expected)
         if i == 2:
             expected = np.array([np.sqrt(0.5), np.sqrt(0.5) * (-1) ** result])
-            np.testing.assert_almost_equal(step.state_vector(), expected)
+            np.testing.assert_almost_equal(step.state_vector(copy=True), expected)
 
 
 @pytest.mark.parametrize('dtype', [np.complex64, np.complex128])
@@ -710,8 +721,8 @@ def _apply_unitary_(self, args: cirq.ApplyUnitaryArgs):
 
     initial_state = np.array([np.sqrt(0.5), np.sqrt(0.5)], dtype=np.complex64)
     result = simulator.simulate(circuit, initial_state=initial_state)
-    np.testing.assert_array_almost_equal(result.state_vector(), initial_state)
-    assert not initial_state is result.state_vector()
+    np.testing.assert_array_almost_equal(result.state_vector(copy=False), initial_state)
+    assert not initial_state is result.state_vector(copy=False)
 
 
 def test_does_not_modify_initial_state():
@@ -735,7 +746,7 @@ def _apply_unitary_(self, args: cirq.ApplyUnitaryArgs):
     result = simulator.simulate(circuit, initial_state=initial_state)
     np.testing.assert_array_almost_equal(np.array([1, 0], dtype=np.complex64), initial_state)
     np.testing.assert_array_almost_equal(
-        result.state_vector(), np.array([0, 1], dtype=np.complex64)
+        result.state_vector(copy=False), np.array([0, 1], dtype=np.complex64)
     )
 
 
@@ -787,7 +798,7 @@ def test_simulates_composite():
     np.testing.assert_allclose(
         c.final_state_vector(ignore_terminal_measurements=False, dtype=np.complex64), expected
     )
-    np.testing.assert_allclose(cirq.Simulator().simulate(c).state_vector(), expected)
+    np.testing.assert_allclose(cirq.Simulator().simulate(c).state_vector(copy=False), expected)
 
 
 def test_simulate_measurement_inversions():
@@ -804,7 +815,7 @@ def test_works_on_pauli_string_phasor():
     a, b = cirq.LineQubit.range(2)
     c = cirq.Circuit(np.exp(0.5j * np.pi * cirq.X(a) * cirq.X(b)))
     sim = cirq.Simulator()
-    result = sim.simulate(c).state_vector()
+    result = sim.simulate(c).state_vector(copy=False)
     np.testing.assert_allclose(result.reshape(4), np.array([0, 0, 0, 1j]), atol=1e-8)
 
 
@@ -812,7 +823,7 @@ def test_works_on_pauli_string():
     a, b = cirq.LineQubit.range(2)
     c = cirq.Circuit(cirq.X(a) * cirq.X(b))
     sim = cirq.Simulator()
-    result = sim.simulate(c).state_vector()
+    result = sim.simulate(c).state_vector(copy=False)
     np.testing.assert_allclose(result.reshape(4), np.array([0, 0, 0, 1]), atol=1e-8)
 
 
@@ -1322,9 +1333,9 @@ def test_final_state_vector_is_not_last_object():
     initial_state = np.array([1, 0], dtype=np.complex64)
     circuit = cirq.Circuit(cirq.wait(q))
     result = sim.simulate(circuit, initial_state=initial_state)
-    assert result.state_vector() is not initial_state
-    assert not np.shares_memory(result.state_vector(), initial_state)
-    np.testing.assert_equal(result.state_vector(), initial_state)
+    assert result.state_vector(copy=False) is not initial_state
+    assert not np.shares_memory(result.state_vector(copy=False), initial_state)
+    np.testing.assert_equal(result.state_vector(copy=False), initial_state)
 
 
 def test_deterministic_gate_noise():
diff --git a/cirq-core/cirq/sim/state_vector.py b/cirq-core/cirq/sim/state_vector.py
index fcd8e62e02e..928722dc0f2 100644
--- a/cirq-core/cirq/sim/state_vector.py
+++ b/cirq-core/cirq/sim/state_vector.py
@@ -57,7 +57,7 @@ def _qid_shape_(self) -> Tuple[int, ...]:
         return self._qid_shape
 
     @abc.abstractmethod
-    def state_vector(self) -> np.ndarray:
+    def state_vector(self, copy: Optional[bool] = False) -> np.ndarray:
         """Return the state vector (wave function).
 
         The vector is returned in the computational basis with these basis
@@ -83,6 +83,12 @@ def state_vector(self) -> np.ndarray:
                 |  6  |   1    |   1    |   0    |
                 |  7  |   1    |   1    |   1    |
 
+        Args:
+            copy: If True, the returned state vector will be a copy of that
+            stored by the object. This is potentially expensive for large
+            state vectors, but prevents mutation of the object state, e.g. for
+            operating on intermediate states of a circuit.
+            Defaults to False.
         """
         raise NotImplementedError()
 
@@ -95,7 +101,9 @@ def dirac_notation(self, decimals: int = 2) -> str:
         Returns:
             A pretty string consisting of a sum of computational basis kets
             and non-zero floats of the specified accuracy."""
-        return qis.dirac_notation(self.state_vector(), decimals, qid_shape=self._qid_shape)
+        return qis.dirac_notation(
+            self.state_vector(copy=False), decimals, qid_shape=self._qid_shape
+        )
 
     def density_matrix_of(self, qubits: List['cirq.Qid'] = None) -> np.ndarray:
         r"""Returns the density matrix of the state.
@@ -132,7 +140,7 @@ def density_matrix_of(self, qubits: List['cirq.Qid'] = None) -> np.ndarray:
                 corresponding to the state.
         """
         return qis.density_matrix_from_state_vector(
-            self.state_vector(),
+            self.state_vector(copy=False),
             [self.qubit_map[q] for q in qubits] if qubits is not None else None,
             qid_shape=self._qid_shape,
         )
@@ -157,7 +165,7 @@ def bloch_vector_of(self, qubit: 'cirq.Qid') -> np.ndarray:
                 corresponding to the state.
         """
         return qis.bloch_vector_from_state_vector(
-            self.state_vector(), self.qubit_map[qubit], qid_shape=self._qid_shape
+            self.state_vector(copy=False), self.qubit_map[qubit], qid_shape=self._qid_shape
         )
 
 
diff --git a/cirq-core/cirq/sim/state_vector_simulator.py b/cirq-core/cirq/sim/state_vector_simulator.py
index daf3dc7a738..44c9aa1a10c 100644
--- a/cirq-core/cirq/sim/state_vector_simulator.py
+++ b/cirq-core/cirq/sim/state_vector_simulator.py
@@ -29,7 +29,7 @@
 
 import numpy as np
 
-from cirq import ops, value, qis
+from cirq import _compat, ops, value, qis
 from cirq._compat import deprecated_class, proper_repr
 from cirq.sim import simulator, state_vector, simulator_base
 
@@ -168,7 +168,7 @@ def final_state_vector(self) -> np.ndarray:
             self._final_state_vector = tensor
         return self._final_state_vector
 
-    def state_vector(self) -> np.ndarray:
+    def state_vector(self, copy: bool = None) -> np.ndarray:
         """Return the state vector at the end of the computation.
 
         The state is returned in the computational basis with these basis
@@ -194,7 +194,13 @@ def state_vector(self) -> np.ndarray:
                 |  6  |   1    |   1    |   0    |
                 |  7  |   1    |   1    |   1    |
         """
-        return self.final_state_vector.copy()
+        if copy is None:
+            _compat._warn_or_error(
+                "Starting in v0.16, state_vector will not copy the state by default. "
+                "Explicitly set copy=True to copy the state."
+            )
+            copy = True
+        return self.final_state_vector.copy() if copy else self.final_state_vector
 
     def _value_equality_values_(self):
         measurements = {k: v.tolist() for k, v in sorted(self.measurements.items())}
diff --git a/cirq-core/cirq/sim/state_vector_simulator_test.py b/cirq-core/cirq/sim/state_vector_simulator_test.py
index 6bff10df026..b98353a544a 100644
--- a/cirq-core/cirq/sim/state_vector_simulator_test.py
+++ b/cirq-core/cirq/sim/state_vector_simulator_test.py
@@ -140,7 +140,23 @@ def test_state_vector_trial_state_vector_is_copy():
     trial_result = cirq.StateVectorTrialResult(
         params=cirq.ParamResolver({}), measurements={}, final_simulator_state=final_simulator_state
     )
-    assert trial_result.state_vector() is not final_simulator_state.target_tensor
+    assert trial_result.state_vector(copy=True) is not final_simulator_state.target_tensor
+
+
+def test_implicit_copy_deprecated():
+    final_state_vector = np.array([0, 1], dtype=np.complex64)
+    qubit_map = {cirq.NamedQubit('a'): 0}
+    final_simulator_state = cirq.StateVectorSimulationState(
+        qubits=list(qubit_map), initial_state=final_state_vector
+    )
+    trial_result = cirq.StateVectorTrialResult(
+        params=cirq.ParamResolver({}), measurements={}, final_simulator_state=final_simulator_state
+    )
+
+    with cirq.testing.assert_deprecated(
+        "state_vector will not copy the state by default", deadline="v0.16"
+    ):
+        _ = trial_result.state_vector()
 
 
 def test_str_big():
diff --git a/cirq-core/cirq/sim/state_vector_test.py b/cirq-core/cirq/sim/state_vector_test.py
index d8b4df663de..4f5f2ea342c 100644
--- a/cirq-core/cirq/sim/state_vector_test.py
+++ b/cirq-core/cirq/sim/state_vector_test.py
@@ -14,6 +14,7 @@
 """Tests for state_vector.py"""
 
 import itertools
+from typing import Optional
 import pytest
 
 import numpy as np
@@ -24,7 +25,7 @@
 
 def test_state_mixin():
     class TestClass(cirq.StateVectorMixin):
-        def state_vector(self) -> np.ndarray:
+        def state_vector(self, copy: Optional[bool] = None) -> np.ndarray:
             return np.array([0, 0, 1, 0])
 
     qubits = cirq.LineQubit.range(2)
@@ -330,7 +331,7 @@ def test_measure_state_empty_state():
 
 
 class BasicStateVector(cirq.StateVectorMixin):
-    def state_vector(self) -> np.ndarray:
+    def state_vector(self, copy: Optional[bool] = None) -> np.ndarray:
         return np.array([0, 1, 0, 0])
 
 
diff --git a/cirq-google/cirq_google/calibration/engine_simulator.py b/cirq-google/cirq_google/calibration/engine_simulator.py
index ef53ddb3edd..620743cd223 100644
--- a/cirq-google/cirq_google/calibration/engine_simulator.py
+++ b/cirq-google/cirq_google/calibration/engine_simulator.py
@@ -372,7 +372,7 @@ def create_from_characterizations_sqrt_iswap(
 
     def final_state_vector(self, program: cirq.Circuit) -> np.ndarray:
         result = self.simulate(program)
-        return result.state_vector()
+        return result.state_vector(copy=False)
 
     def get_calibrations(
         self, requests: Sequence[PhasedFSimCalibrationRequest]
diff --git a/docs/simulation.ipynb b/docs/simulation.ipynb
index 691f125f0f7..843d8ee5fad 100644
--- a/docs/simulation.ipynb
+++ b/docs/simulation.ipynb
@@ -518,7 +518,7 @@
     "circuit = cirq.Circuit()\n",
     "circuit.append(basic_circuit())    \n",
     "for i, step in enumerate(simulator.simulate_moment_steps(circuit)):\n",
-    "    print('state at step %d: %s' % (i, np.around(step.state_vector(), 3)))"
+    "    print('state at step %d: %s' % (i, np.around(step.state_vector(copy=True), 3)))"
    ]
   },
   {
diff --git a/docs/tutorials/educators/intro.ipynb b/docs/tutorials/educators/intro.ipynb
index ca371aa4f7e..12634b3a2d0 100644
--- a/docs/tutorials/educators/intro.ipynb
+++ b/docs/tutorials/educators/intro.ipynb
@@ -1811,7 +1811,7 @@
     "\n",
     "# Step through the simulation results.\n",
     "for step in simulator.simulate_moment_steps(circuit):\n",
-    "    prob = np.abs(step.state_vector()) ** 2\n",
+    "    prob = np.abs(step.state_vector(copy=True)) ** 2\n",
     "    probs.append(prob[0])\n",
     "\n",
     "# Plot the probability of the ground state at each simulation step.\n",
diff --git a/docs/tutorials/fourier_checking.ipynb b/docs/tutorials/fourier_checking.ipynb
index 08c06f419e8..5b8887f9768 100644
--- a/docs/tutorials/fourier_checking.ipynb
+++ b/docs/tutorials/fourier_checking.ipynb
@@ -743,7 +743,7 @@
     "for step in s.simulate_moment_steps(circuit):\n",
     "    print(step.dirac_notation())\n",
     "    print(\"|0> state probability to observe: \",\n",
-    "          np.abs(step.state_vector()[0])**2)"
+    "          np.abs(step.state_vector(copy=True)[0])**2)"
    ]
   },
   {
diff --git a/examples/stabilizer_code.ipynb b/examples/stabilizer_code.ipynb
index 63edf1492f6..d2e99448681 100644
--- a/examples/stabilizer_code.ipynb
+++ b/examples/stabilizer_code.ipynb
@@ -101,7 +101,7 @@
    "source": [
     "for initial_state in [0, 1]:\n",
     "    results = cirq.Simulator().simulate(circuit, qubit_order=qubits, initial_state=initial_state)\n",
-    "    print('%d becomes %s' % (initial_state, cirq.dirac_notation(results.state_vector())))\n",
+    "    print('%d becomes %s' % (initial_state, cirq.dirac_notation(results.state_vector(copy=False))))\n",
     "    \n",
     "    "
    ]
@@ -160,7 +160,7 @@
    ],
    "source": [
     "results = cirq.Simulator().simulate(circuit, qubit_order=(qubits+ancillas), initial_state=0)\n",
-    "print('%s' % (cirq.dirac_notation(results.state_vector())))"
+    "print('%s' % (cirq.dirac_notation(results.state_vector(copy=False))))"
    ]
   },
   {
@@ -216,7 +216,7 @@
     }
    ],
    "source": [
-    "decoded = code.decode(qubits, ancillas, results.state_vector())\n",
+    "decoded = code.decode(qubits, ancillas, results.state_vector(copy=False))\n",
     "print('decoded=%s' % (decoded))"
    ]
   },
@@ -256,7 +256,7 @@
     "\n",
     "for initial_state in [0, 1]:\n",
     "    results = cirq.Simulator().simulate(circuit, qubit_order=qubits, initial_state=initial_state)\n",
-    "    print('%d becomes:\\n %s\\n' % (initial_state, cirq.dirac_notation(4.0 * results.state_vector())))"
+    "    print('%d becomes:\\n %s\\n' % (initial_state, cirq.dirac_notation(4.0 * results.state_vector(copy=False))))"
    ]
   }
  ],
diff --git a/examples/stabilizer_code_test.py b/examples/stabilizer_code_test.py
index ddd28e2b71a..8a8615c73a8 100644
--- a/examples/stabilizer_code_test.py
+++ b/examples/stabilizer_code_test.py
@@ -37,7 +37,7 @@ def encode_corrupt_correct(
         .simulate(
             circuit, qubit_order=(qubits + ancillas), initial_state=(input_val * 2 ** len(ancillas))
         )
-        .state_vector()
+        .state_vector(copy=False)
     )
 
     decoded = code.decode(qubits, ancillas, state_vector)